home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / hayat.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  113.5 KB  |  3,181 lines

  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;   **************************************************************
  9. ;;;   ***** HAYAT ******* Finite Power Series Routines *************
  10. ;;;   **************************************************************
  11. ;;;   ** (c) Copyright 1982 Massachusetts Institute of Technology **
  12. ;;;   ****** This is a read-only file! (All writes reserved) *******
  13. ;;;   **************************************************************
  14.  
  15. (in-package "MAXIMA")
  16.  
  17. ;;;        TOP LEVEL STRUCTURE
  18.  
  19. ;;;    Power series have the following format when seen outside the power
  20. ;;; series package:
  21. ;;; 
  22. ;;;    ((MRAT SIMP <varlist> <genvar> <tlist> trunc) <poly-form>)
  23. ;;; 
  24. ;;; This is the form of the output of the expressions, to
  25. ;;; be displayed they are RATDISREPed and passed to DISPLA.
  26.  
  27. ;;; The <poly-forms> consist of a header and list of exponent-coefficient
  28. ;;; pairs as shown below.  The PS is used to distinguish power series
  29. ;;; from their coefficients which have a similar representation.
  30. ;;; 
  31. ;;;   (PS (<var> . <ord-num>) (<trunc-lvl>)
  32. ;;;      (<exponent> . <coeff>) (<exponent> . <coeff>) . . .)
  33. ;;; 
  34. ;;; The <var> component of the power series is a gensym which represents the
  35. ;;; kernel of the power series.  If the package is called with the arguments:
  36. ;;; Taylor(<expr>, x, a, n)  then the kernel will be (x - a).
  37. ;;; The <ord-num> is a relative ordering for the various kernels in a
  38. ;;; multivariate expansion.  
  39. ;;; <trunc-lvl> is the highest degree of the variable <var> which is retained
  40. ;;; in the current power series.
  41. ;;; The terms in the list of exponent-coefficient pairs are ordered by
  42. ;;; increasing degree.
  43.  
  44. ;;; Problem: fix expansion of logs so that taylor(log(1+exp(-1/x)),x,0,3)
  45. ;;; works. Done.
  46. ;;; 
  47. ;;; Problem: taylor(log(1+exp(-1/x)),x,0,5) loses because, while
  48. ;;; taylor_simplify_recurse'ing exp(-3/x) get trunc level = -3. FIxed.
  49. ;;; 
  50. ;;; Problem: Need to fix things so that asymptotic kernels aren't put onto
  51. ;;; tvars via tlist merge etc. in taylor1. Done.
  52. ;;; 
  53. ;;; Problem: get-series returns 0 for taylor(log(1+exp(-1/x)),x,0,5) and
  54. ;;; need to make log(exp(1/x)) -> 1/x. Fixed.
  55. ;;; 
  56. ;;; Problem: Fix psexpt-fn so that it doesn't lose via the invert-var
  57. ;;; scheme, e.g. try taylor(exp(exp(-1/x)+x),x,0,5). Note that if it did
  58. ;;; just that scheme. Done.
  59. ;;; 
  60. ;;; Problem: fix adjoin-tvar so that the new tvars are ordered correctly
  61. ;;; according to their strength. This is necessary in order to read the limit
  62. ;;; directly from the leading term. E.g. see the misordered:
  63. ;;; taylor(subst(1/x,x,part(screw2,1)),x,0,2) from ALJABR;SCREW2 LIMIT.
  64. ;;; Note that the answer given for this appear to be incorrect when the
  65. ;;; truncation on x is < 4. Is this due to the misordering?
  66. ;;; Also taylor(screwa,x,0,4)+taylor(screwb,x,0,8) doesn't agree with
  67. ;;; taylor(screw,x,0,8) where it should (here screwa = part(screw,1),
  68. ;;; screwb = part(screw, 2); is this a truncation problem on the
  69. ;;; gvar exp(1/x)?).
  70. ;;; 
  71. ;;; Problem: new gvars have to be intro'd for logs just as for exp's instead
  72. ;;; of treating them like constants as currently done. For example,
  73. ;;; taylor(log(1+1/log(x)),x,0,2) currently doesn't expand. Done.
  74. ;;; 
  75. ;;; Problem: The display routines need pieces of the taylor environment
  76. ;;; (tvar-limits, tvars, tlist, etc.) to figure out how to order terms.
  77. ;;; This means we'll probably have to store some of this on the local tlist.
  78. ;;; When this is done the commented out code in psdisrep and psdisrep2 can
  79. ;;; be re-inserted. Psdisrep2expand will also need to be modified.
  80. ;;; I just fixed srdisrep to get the local env it needs; psdisrep2expand
  81. ;;; still needs to be updated. The display order problem is still around 
  82. ;;; however: try taylor(exp(exp(-1/x)+x),x,0,3). After more investigation,
  83. ;;; it seems that the term reversal always occurs for ps's that are coeff's
  84. ;;; of terms whose expt is < 0. Probably the psdisrep routines should reverse
  85. ;;; these terms to account for this (the bug is somewhere in the DISPLA
  86. ;;; routines, possible DIM-MPLUS).
  87. ;;;
  88. ;;; Problem: Since gvar's like exp(-1/x) can't be put on genvar, they have
  89. ;;; to be saved somewhere locally and restored by everyone who needs to setup
  90. ;;; disrep info etc. Done for re-taylor.
  91. ;;;
  92. ;;; Problem: All new code needs to be checked to ensure it does the correct
  93. ;;; thing when the expansion point is infinite (e.g. see the code in
  94. ;;; TSEXPT-RED which handles this case).
  95. ;;;
  96. ;;; Perhaps the code for exp's and log's which pushes trunc degrees
  97. ;;; can be done by first computing exp(c0) or log(c0) first and see
  98. ;;; how much to push by looking at this series. Done for exp in tsexpt-red.
  99. ;;;
  100. ;;; Problems: taylor(part(screwa,2)-2/x,x,0,1) shouldn't be exact.
  101. ;;; taylor(screwa,x,0,-2) misses the degree -2 term. This part is now fixed.
  102. ;;;
  103. ;;; Tvar-limits should be stored locally so that psdisrep need not recompute
  104. ;;; each gvar limit when disrepping.
  105.  
  106. (macsyma-module hayat)
  107.  
  108. (defmvar tlist nil)
  109. (defvar *within-srf?* nil)
  110.  
  111. (load-macsyma-macros mhayat rzmac ratmac)
  112.  
  113. ;        (comment Subtitle Special Stuff for Compiling)
  114.  
  115. (declare-top(special vlist        ;
  116.           varlist    ;List of all the variables occuring in a power
  117.                 ;series, the power series variables at the end
  118.           genvar    ;The list of gensyms corresponding to varlist
  119.           modulus       ;
  120.           *a*        ;Temporary special
  121.           sub-exprs    ;
  122.           silent-taylor-flag    ;If true indicates that errors will be 
  123.                 ;returned via a throw to TAY-ERR
  124.           tlist        ;An association list which contains the
  125.                 ;relevant information for the expansion which
  126.                 ;is passed in at toplevel invocation.  
  127.           $float    ;Indicates whether to convert rational numbers
  128.                 ;to floating point numbers.
  129.           $keepfloat    ;When true retains floatin point numbers 
  130.                 ;internal to Taylor.
  131.           $radexpand    ;
  132.           log-1        ;What log(-1) should be log(-1) or pi*i.
  133.           log%i        ;Similarly for log(i)
  134.           exact-poly    ;Inicates whether polynomials are to be
  135.                         ;considered exact or not.  True within SRF,
  136.                 ;false within TAYLOR.
  137.           ngps        ;
  138.           num-syms    ;
  139.           loc-gensym    ;
  140.           syms        ;
  141.           tvars        ;
  142.           pssyms    ;
  143.           half%pi    ;Has pi/2 to save space.
  144.           const-funs    ;
  145.           const-exp-funs;
  146.           tay-const-expand    ;For rediculousness like csch(log(x))
  147.           $exponentialize       ;which we do by exponentiation.
  148.           tay-pole-expand;
  149.           trigdisp    ;
  150.           last-exp    ;last-expression through taylor2
  151.           $taylordepth    ;
  152.           $ratexpand    ;
  153.           genpairs    ;List of dotted pairs 
  154.           ps-bmt-disrep    ;
  155.           ivars        ;Pairlist if gensym and disreped version
  156.           key-vars    ;Pairlist of gensym and key var (for searching
  157.                 ;TLIST)
  158.           $algebraic    ;
  159.           *psacirc    ;
  160.           *pscirc    ;
  161.           full-log    ;
  162.           $logarc    ;
  163.           trunclist    ;
  164.           *within-srf?* ;flag for in srf
  165.           mainvar-datum ;
  166.           least_term?    ; If non-null then the addition routines
  167.                 ; are adding or subtracting coeff's of the
  168.                 ; least term of a sum so they should do
  169.                 ; zero checking on it if it is desired.
  170.           taylor_simplifier ; This is set by taylor1 to the object
  171.                     ; which will be funcalled whenever
  172.                     ; coefficient simplification is desired.
  173.           zerolist      ; A list of constant expressions which have
  174.                 ; been verified to be zero by a call to
  175.                 ; $TAYLOR_SIMPLIFIER in taylor2. It is used to
  176.                 ; suppress the message that TAYLOR is assumming
  177.                 ; an expression to be zero.
  178.               $zerobern $simp 0p-funord lexp-non0)    ;
  179.      (*expr lcm)
  180.      (muzzled t))    ;Don't want to see closed compilation notes.
  181.  
  182. (defmvar $PSEXPAND ()
  183.  "When TRUE extended rational function expressions will be displayed fully
  184.   expanded. (RATEXPAND will also cause this.) If FALSE, multivariate
  185.   expressions will be displayed just as in the rational function package.
  186.   If PSEXPAND:MULTI, then terms with the same total degree in the variables
  187.   are grouped together.")
  188.  
  189. (defmvar $MAXTAYORDER T
  190.  "When true TAYLOR retains as many terms as are certain to be correct
  191.   during power series arithmetic. Otherwise, truncation is controlled
  192.   by the arguments specified to TAYLOR.")
  193.  
  194. (defmvar $TAYLOR_TRUNCATE_POLYNOMIALS T
  195.  "When FALSE polynomials input to TAYLOR are considered to have infinite
  196.   precison; otherwise (the default) they are truncated based upon the input
  197.   truncation levels.")
  198.  
  199. (defmvar $TAYLOR_LOGEXPAND ()
  200.  "Unless FALSE log's of products will be expanded fully in TAYLOR (the default)
  201.   to avoid identically-zero constant terms which involve log's. When FALSE,
  202.   only expansions necessary to produce a formal series will be executed.")
  203.  
  204. ;Note!  The value of this must be a symbol, because it is checked with
  205. ; FBOUNDP.
  206. (defmvar $TAYLOR_SIMPLIFIER 'SIMPLIFY
  207.  "A function of one argument which TAYLOR uses to simplify coefficients
  208.   of power series.")
  209.  
  210. ;        (Comment Subtitle General Macsyma Free Predicates)
  211.  
  212. (defun zfree (e x)
  213.     (cond ((equal e x) () )
  214.       ((atom e) 'T)
  215.       ((eq (caar e) 'MRAT)
  216.        (null (zl-MEMBER x (cdr ($listofvars e)))))
  217.       ('T (do ((l (cdr e) (cdr l))) ((null l) 'T)
  218.          (or (zfree (car l) x) (return () ))))))
  219.  
  220. (defun mfree (exp varl)
  221.   (declare (special DUMMY-VARIABLE-OPERATORS))
  222.    (cond ((atom exp) (not (memq exp varl)))
  223.      ((eq (caar exp) 'MRAT)
  224.       (do ((l (mrat-varlist exp) (cdr l)))
  225.           ((null l) 'T)
  226.          (unless (mfree (car l) varl) (return () ))))
  227.      ((or (memq (caar exp) DUMMY-VARIABLE-OPERATORS)
  228.           (memq 'array (cdar exp)))
  229.       (do ((vars varl (cdr vars)))
  230.           ((null vars) 'T)
  231.          (unless (freeof (car vars) exp) (return () ))))
  232.      ('T (and (mfree (caar exp) varl) (mfreel (cdr exp) varl)))))
  233.  
  234. (defun mfreel (l varl)
  235.   (or (null l) (and (mfree (car l) varl) (mfreel (cdr l) varl))))
  236.  
  237. ;(comment Subtitle Coefficient Arithmetic)
  238.  
  239. (defun rcexpt (x y)
  240.        (cond ((equal x (rcone)) (rcone))
  241.          ((rczerop y) (rcone))
  242.          ((and (equal (cdr y) 1) (fixnump (car y)))
  243.           (ratexpt x (car y)))
  244.          ((and $radexpand (numberp (car y)) (numberp (cdr y)))
  245.           (if (floatp (car y))
  246.           (setq y (MAXIMA-RATIONALIZE (*quo (car y) (cdr y)))))
  247.           (ratexpt (rcquo (rcexpt1 (car x) (cdr y))
  248.                   (rcexpt1 (cdr x) (cdr y)))
  249.                (car y)))
  250.          (t (let ($keepfloat)
  251.              (prep1 (m^ (rcdisrep x) (rcdisrep y)))))))
  252.  
  253. (defun rcexpt1 (p n)
  254.    (cond ((equal p 1) (rcone))
  255.      ((pcoefp p) (prep1 (m^ (pdis p) (*red 1 n))))
  256.      ;; psfr does a square-free decom on p yielding (p1 e1 p2 e2 ... pn en)
  257.      ;; where p = p1^e1 p2^e2 ... pn^en, the pi being square-free
  258.      (T (do ((l (psqfr p) (cddr l))
  259.          (ans (rcone)))
  260.         ((null l) ans)
  261.            (if (not (equal (remainder (cadr l) n) 0))
  262.            (setq ans (rctimes ans (prep1 (m^ (pdis (car l))
  263.                              (*red (cadr l) n)))))
  264.           ;; If pi<0, n=2m and n|ei then ei=2e and
  265.           ;;        (pi^ei)^(1/(2m)) = (-pi)^(e/m)
  266.           (when (and (evenp n) (eq ($sign (pdis (car l))) '$NEG))
  267.              (rplaca l (pminus (car l))))
  268.           (setq ans (rctimes ans (ratexpt (cons (car l) 1)
  269.                            (// (cadr l) n)))))))))
  270.  
  271. (defun rccoefp (e)        ;a sure check, but expensive
  272.        (and (null (atom e))
  273.         (or (atom (car e))
  274.         (memq (caar e) genvar))
  275.         (or (atom (cdr e))
  276.         (memq (cadr e) genvar))))
  277.  
  278. ;        (Comment Subtitle Exponent arithmetic)
  279.  
  280. (defun ezerop (x) (and (not (infp x)) (signp e (car x))))
  281.  
  282. (defun e+ (x y)
  283.     (cond ((or (infp x) (infp y)) (inf))
  284.       ((and (equal (cdr x) 1) (equal (cdr y) 1))
  285.        (cons (plus (car x) (car y)) 1))
  286.       (t (ereduce (plus (times (car x) (cdr y)) (times (cdr x) (car y)))
  287.               (times (cdr x) (cdr y))))))
  288.  
  289. (defun ediff (x y)
  290.     (cond ((infp x) (inf))
  291.       ((and (equal (cdr x) 1) (equal (cdr y) 1))
  292.        (cons (*dif (car x) (car y)) 1))
  293.       (t (ereduce (*dif (times (car x) (cdr y)) (times (cdr x) (car y)))
  294.               (times (cdr x) (cdr y))))))
  295.  
  296. (defun emin (x y)
  297.     (cond ((infp x) y)
  298.       ((infp y) x)
  299.       ((equal (cdr x) (cdr y)) (cons (min (car x) (car y)) (cdr x)))
  300.       ((lessp (times (car x) (cdr y)) (times (cdr x) (car y))) x)
  301.       (t y)))
  302.  
  303. (defun emax (x y)
  304.     (cond ((or (infp x) (infp y)) (inf))
  305.       ((equal (cdr x) (cdr y)) (cons (max (car x) (car y)) (cdr x)))
  306.       ((greaterp (times (car x) (cdr y)) (times (cdr x) (car y))) x)
  307.       (t y)))
  308.  
  309. (defun e* (x y)
  310.     (cond ((or (infp x) (infp y)) (inf))
  311.       ((and (equal (cdr x) 1) (equal (cdr y) 1))
  312.        (cons (times (car x) (car y)) 1))
  313.       (t (ereduce (times (car x) (car y)) (times (cdr x) (cdr y))))))
  314.  
  315. (defun erecip (e)
  316.        (if (minusp (car e))
  317.        (cons (minus (cdr e)) (minus (car e)))
  318.        (cons (cdr e) (car e))))
  319.  
  320. (defun equo (x y)
  321.        (cond ((infp x) (inf))
  322.          ((infp y) (rczero))
  323.          (t (ereduce (times (car x) (cdr y))
  324.              (times (cdr x) (car y))))))
  325.  
  326. (defun e1+ (x)
  327.     (cond ((infp x) (inf))
  328.       ((= (cdr x) 1) (cons (add1 (car x)) 1))
  329.       (t (cons (plus (cdr x) (car x)) (cdr x)))))
  330.  
  331. (defun e1- (x)
  332.     (cond ((infp x) (inf))
  333.       ((equal (cdr x) 1) (cons (sub1 (car x)) 1))
  334.       (t (cons (*dif (car x) (cdr x)) (cdr x)))))
  335.  
  336. (defun e> (x y)
  337.     (cond ((infp x) t)
  338.       ((infp y) ())
  339.       ((equal (cdr x) (cdr y)) (greaterp (car x) (car y)))
  340.       (t (greaterp (times (car x) (cdr y)) (times (car y) (cdr x))))))
  341.  
  342. (defun e= (e1 e2)
  343.       (cond ((eq e1 e2) t)
  344.         ((or (null e1) (null e2)) ())
  345.         (t (and (equal (car e1) (car e2))
  346.             (equal (cdr e1) (cdr e2))))))
  347.  
  348. (defun ereduce (n d)
  349.        (if (signp l d) (setq d (minus d) n (minus n)))
  350.        (if (zerop n) (rczero)
  351.        (let ((gcd (gcd n d)))
  352.         (cons (*quo n gcd) (*quo d gcd)))))
  353.  
  354. (defun egcd (x y)
  355.        (let ((xn (abs (car x))) (xd (cdr x))
  356.          (yn (abs (car y))) (yd (cdr y)))
  357.         (cons (gcd xn yn) (times xd (*quo yd (gcd xd yd))))))
  358.  
  359. ;        (Comment Subtitle polynomial arithmetic)
  360.  
  361. (declare-top(special vars))
  362.  
  363. (defun ord-vector (p)
  364.    (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum)))
  365.                tlist)))
  366.      (declare (special vars))
  367.       (ifn (cdr vars) (ncons (ps-le* p))
  368.      (ord-vect1 p)
  369.      (mapcar #'(lambda (x) (or (cdr x) (rczero)))
  370.          vars))))
  371.  
  372. (defun ord-vect1 (p)
  373.   (declare (special vars))
  374.   (unless (pscoefp p)
  375.      (let ((data (assq (gvar p) vars))
  376.        (le (ps-le p)))
  377.     (rplacd data (ifn (cdr data) le (emin (cdr data) le)))
  378.     (MAPL #'(lambda (l) (ord-vect1 (lc l)))
  379.          (terms p)))))
  380.  
  381. (defun trunc-vector (p min?)
  382.    (let ((vars (mapcar #'(lambda (datum) (list (int-gvar datum)))
  383.                tlist)))
  384.      (declare (special vars))
  385.       (if (null (cdr vars)) (ncons (if (psp p) (trunc-lvl p) () ))
  386.      (trunc-vect1 p min?)
  387.      (mapcar 'cdr vars))))
  388.  
  389. (defun trunc-vect1 (p min?)
  390.   (declare (special vars))
  391.    (unless (pscoefp p)
  392.       (let ((data (assq (gvar p) vars))
  393.         (trunc (trunc-lvl p)))
  394.      (when trunc
  395.         (rplacd data (if (null (cdr data)) trunc
  396.                (if min? (emin (cdr data) trunc)
  397.                   (emax (cdr data) trunc))))))
  398.       (dolist (term (terms p)) (trunc-vect1 (c term) min?))))
  399.  
  400. (declare-top(unspecial vars))
  401.  
  402. (defun psplus (x y)
  403.    (cond ((pscoefp x)
  404.       (cond ((pscoefp y) (rcplus x y))
  405.         ((rczerop x) y)
  406.         (T (pscplus x y))))
  407.      ((pscoefp y) (if (rczerop y) x (pscplus y x)))
  408.      ((eqgvar (gvar-o x) (gvar-o y)) (psplus1 x y))
  409.      ((pointerp (gvar-o x) (gvar-o y)) (pscplus y x))
  410.      (T (pscplus x y))))
  411.  
  412. (defun rcplus! (x y)
  413.    (if (not (and least_term? taylor_simplifier)) (rcplus x y)
  414.       (prep1 (funcall taylor_simplifier (m+ (rcdisrep x) (rcdisrep y))))))
  415.  
  416. (defun psdiff (x y)
  417.    (cond ((pscoefp x) (cond ((pscoefp y) (rcdiff x y))
  418.                 ((rczerop x) (pstimes (rcmone) y))
  419.                 (T (pscdiff x y () ))))
  420.      ((pscoefp y) (if (rczerop y) x (pscdiff y x T)))
  421.      ((eqgvar (gvar-o x) (gvar-o y)) (psdiff1 x y))
  422.      ((pointerp (gvar-o x) (gvar-o y)) (pscdiff y x T))
  423.      (T (pscdiff x y () ))))
  424.  
  425. (defun rcdiff! (x y)
  426.    (if (not (and least_term? taylor_simplifier)) (rcdiff x y)
  427.       (prep1 (funcall taylor_simplifier (m- (rcdisrep x) (rcdisrep y))))))
  428.  
  429. (defun psplus1 (x y)
  430.    (let ((ans (cons () () )))
  431.       (psplus2 (gvar-o x) (emin (trunc-lvl x) (trunc-lvl y))
  432.            (cons 0 (terms x)) (cons 0 (terms y)) ans ans)))
  433.  
  434. (defun pscplus (c p)
  435.    (if (e> (rczero) (trunc-lvl p)) p
  436.       (pscheck (gvar-o p) (poly-data p) (pscplus1 c (terms p)))))
  437.  
  438. (defun pscdiff (c p fl)
  439.    (if (e> (rczero) (trunc-lvl p))
  440.        (if fl p (psminus p))
  441.       (pscheck (gvar-o p) (poly-data p)
  442.            (ifn fl (pscplus1 c (psminus-terms (terms p)))
  443.           (pscplus1 (psminus c) (terms p))))))
  444.  
  445. (defun strip-zeroes (terms ps?)
  446.    (cond ((or (null terms) (null taylor_simplifier)) terms)
  447.      ((null ps?)
  448.       (do ((terms terms (n-term terms)))
  449.           ((null terms) () )
  450.          (change-lc terms (strip-zeroes (lc terms) 'T))
  451.          (unless (rczerop (lc terms)) (return terms))))
  452.      ((pscoefp terms) 
  453.       (if (null taylor_simplifier) terms
  454.          (let ((exp (rcdisrep terms)))
  455.         ;; If a pscoeff is not free of tvars then the ps is a
  456.         ;; multivar series and we can't handle a recursive
  457.         ;; call to taylor (as opposed to a call to prep1, as below)
  458.         ;; because this would be circuler (e.g. try
  459.         ;; taylor(x/ (x^2+1),[x],%i,-1) ). Besides, in this case
  460.         ;; the pscoeff contains a tvar hence should not be 0.
  461.         (if (not (mfree exp tvars)) terms
  462.            (prep1 (funcall taylor_simplifier exp))))))
  463.      (T (pscheck (gvar-o terms) (poly-data terms)
  464.              (strip-zeroes (terms terms) () )))))
  465.  
  466. (defun pscplus1 (c l)
  467.    (cond ((null l) (list (term (rczero) c)))
  468.      ((rczerop (le l)) (setq c (psplus c (lc l)))
  469.       (if (rczerop c) (strip-zeroes (n-term l) () )
  470.          (cons (term (rczero) c) (n-term l))))
  471.      ((e> (le l) (rczero)) (cons (term (rczero) c) l))
  472.      (T (cons (lt l) (let ((least_term?)) (pscplus1 c (n-term l)))))))
  473.  
  474. ;;; Both here and in psdiff2 xx and yy point one before where one
  475. ;;; might think they should point so that extensions will be retained.
  476.  
  477. (defun psplus2 (varh trunc xx yy ans a)
  478.   (prog (c)
  479.    a    (cond ((mono-term? xx)
  480.            (if (mono-term? yy) (go end) (go null)))
  481.           ((mono-term? yy) (setq yy xx) (go null)))
  482.         (cond ((equal (le (n-term xx)) (le (n-term yy)))
  483.            (setq xx (n-term xx) yy (n-term yy))
  484.            (setq c (let ((least_term? (null (n-term ans))))
  485.               (psplus (lc xx) (lc yy))))
  486.            (if (rczerop c) (go a) (add-term a (le xx) c)))
  487.           ((e> (le (n-term xx)) (le (n-term yy)))
  488.            (setq yy (n-term yy))
  489.            (add-term a (lt yy)))
  490.           (T (setq xx (n-term xx))
  491.          (add-term a (lt xx))))
  492.     (setq a (n-term a))
  493.     (go a)
  494.    null (if (or (mono-term? yy) (e> (le (n-term yy)) trunc))
  495.         (go end)
  496.         (setq yy (n-term yy))
  497.         (add-term-&-pop a (lt yy))
  498.         (go null))
  499.    end  (return (pscheck varh (list trunc) (cdr ans)))))
  500.  
  501. (defun psdiff1 (x y)
  502.    (let ((ans (cons () () )))
  503.       (psdiff2 (gvar-o x) (emin (trunc-lvl x) (trunc-lvl y))
  504.            (cons 0 (terms x)) (cons 0 (terms y)) ans ans)))
  505.  
  506. (defun psdiff2 (varh trunc xx yy ans a)
  507.   (prog (c)
  508.    a    (cond ((mono-term? xx)
  509.            (if (mono-term? yy) (go end)
  510.            (setq yy (cons 0 (mapcar
  511.                      #'(lambda (q)
  512.                       (term (e q) (psminus (c q))))
  513.                      (cdr yy))))
  514.            (go null)))
  515.           ((mono-term? yy)
  516.            (setq yy xx) (go null)))
  517.         (cond ((equal (le (n-term xx)) (le (n-term yy)))
  518.            (setq xx (n-term xx) yy (n-term yy))
  519.            (setq c (let ((least_term? (null (n-term ans))))
  520.               (psdiff (lc xx) (lc yy))))
  521.            (if (rczerop c) (go a)
  522.            (add-term a (le xx) c)))
  523.           ((e> (le (n-term xx)) (le (n-term yy)))
  524.            (setq yy (n-term yy))
  525.            (add-term a (le yy) (psminus (lc yy))))
  526.           (T (setq xx (n-term xx))
  527.          (add-term a (lt xx))))
  528.     (setq a (n-term a))
  529.     (go a)
  530.    null (if (or (mono-term? yy) (e> (le (n-term yy)) trunc))
  531.         (go end)
  532.         (setq yy (n-term yy))
  533.         (add-term-&-pop a (le yy) (lc yy))
  534.         (go null))
  535.    end    (return (pscheck varh (list trunc) (cdr ans)))))
  536.  
  537. (defun psminus (x)
  538.    (if (psp x) (make-ps x (psminus-terms (terms x)))
  539.       (rcminus x)))
  540.  
  541. (defun psminus-terms (terms)
  542.    (let ((ans (cons () () )))
  543.       (do ((q terms (n-term q))
  544.        (a ans (cdr a)))
  545.       ((null q) (cdr ans))
  546.      (add-term a (le q) (psminus (lc q))))))
  547.  
  548. (defun pscheck (a b terms)
  549.    (cond ((null terms) (rczero))
  550.      ((and (mono-term? terms) (rczerop (le terms)))
  551.       (lc terms))
  552.      (T (make-ps a b terms))))
  553.  
  554. (defun pstrim-terms (terms e)
  555.    (do () (())
  556.       (cond ((null terms) (return () ))
  557.         ((null (e> e (le terms))) (return terms))
  558.         (T (setq terms (n-term terms))))))
  559.  
  560. (defun psterm (terms e)
  561.    (psterm1 (pstrim-terms terms e) e))
  562.  
  563. (defun psterm1 (l e) 
  564.    (cond ((null l) (rczero))
  565.      ((e= (le l) e) (lc l))
  566.      (T (rczero))))
  567.  
  568. (defun pscoeff1 (a b c)        ;a is an mrat!!!
  569.    (let ((tlist (mrat-tlist a)))
  570.       (cons (nconc (list 'MRAT 'SIMP (mrat-varlist a) (mrat-genvar a))
  571.            (do ((l (mrat-tlist a) (cdr l))
  572.             (ans () (cons (car l) ans)))
  573.                ((null l) ans)
  574.               (when (alike1 (caar l) b)
  575.              (return
  576.               (and (or ans (cdr l))
  577.                    (list (nreconc ans (cdr l)) 'TRUNC))))))
  578.         (pscoef (mrat-ps a) (int-gvar (get-datum b)) (prep1 c)))))
  579.  
  580. (defun pscoef (a b c)
  581.    (cond ((pscoefp a) (if (rczerop c) a (rczero)))
  582.      ((eq b (gvar a)) (psterm (terms a) c))
  583.      (T (do ((gvar-o (gvar-o a))
  584.          (poly-data (poly-data a))
  585.          (ans (rczero))
  586.          (terms (terms a) (n-term terms))
  587.          (temp))
  588.         ((null terms) ans)
  589.            (unless (rczerop (setq temp (pscoef (lc terms) b c)))
  590.           (setq ans (psplus ans
  591.                     (make-ps gvar-o poly-data
  592.                          (ncons (term (le terms)
  593.                               temp))))))))))
  594.  
  595. (defun psdisextend (p)
  596.   (ifn (psp p) p
  597.        (make-ps p
  598.         (mapcar #'(lambda (q) (cons (car q)
  599.                         (psdisextend (cdr q))))
  600.             (terms p)))))
  601.  
  602. (defun psfloat (p)
  603.    (if (psp p) (psfloat1 p (trunc-lvl p) (terms p) (ncons 0))
  604.       (rctimes (rcfone) p)))
  605.  
  606. (defun psfloat1 (p trunc l ans)
  607.    (do (($float 'T)
  608.     (a (last ans) (n-term a)))
  609.        ((or (null l) (e> (le l) trunc))
  610.     (pscheck (gvar-o p) (poly-data p) (cdr ans)))
  611.       (add-term a (le l) (psfloat (lc l)))
  612.       (setq l (n-term l))))
  613.  
  614. (defun pstrunc (p)
  615.   (pstrunc1 p (mapcar #'(lambda (q) (cons (int-gvar q) (current-trunc q)))
  616.               tlist)))
  617.  
  618. (defun pstrunc1 (p trlist)
  619.    (ifn (psp p) p
  620.       (let ((trnc (cdr (assq (gvar p) trlist))) (trunc-ps) (a () ))
  621.      (do ((l (terms p) (n-term l)))
  622.          ((null l) (pscheck (gvar-o p) (ncons (trunc-lvl p)) (nreverse a)))
  623.         (when (e> (le l) trnc)
  624.            (return (pscheck (gvar-o p) (ncons trnc) (nreverse a))))
  625.         (unless (rczerop (setq trunc-ps (pstrunc1 (lc l) trlist)))
  626.            (push (term (le l) trunc-ps) a))))))
  627.  
  628. (defun pstimes (x y)
  629.    (cond ((or (rczerop x) (rczerop y)) (rczero))
  630.      ((pscoefp x) (cond ((pscoefp y) (rctimes x y))
  631.                 ((equal x (rcone)) y)
  632.                 (T (psctimes* x y))))
  633.      ((pscoefp y) (if (equal y (rcone)) x (psctimes* y x)))
  634.      ((eqgvar (gvar-o x) (gvar-o y)) (pstimes*1 x y))
  635.      ((pointerp (gvar-o x) (gvar-o y)) (psctimes* y x))
  636.      (T (psctimes* x y))))
  637.  
  638. (defun psctimes* (c p)
  639.   (make-ps p (maplist #'(lambda (l)
  640.                (term (le l) (pstimes c (lc l))))
  641.               (terms p))))
  642.  
  643. (defun pstimes*1 (xa ya)
  644.    (let ((ans (cons () () ))
  645.      (trunc (let ((lex (ps-le xa)) (ley (ps-le ya)))
  646.            (e+ (emin (e- (trunc-lvl xa) lex) (e- (trunc-lvl ya) ley))
  647.                (e+ lex ley)))))
  648.       (unless $maxtayorder
  649.      (setq trunc (emin trunc (t-o-var (gvar xa)))))
  650.       (pstimes*2 xa ya trunc ans)))
  651.  
  652. (defun pstimes*2 (xa ya trunc ans)
  653.    (prog (a c e x y yy)
  654.      (setq x (terms xa) y (setq yy (terms ya)) a ans)
  655.     a     (cond ((or (null y) (e> (setq e (e+ (le x) (le y))) trunc))
  656.         (go b))
  657.            ((not (rczerop (setq c (pstimes (lc x) (lc y)))))
  658.         (add-term-&-pop a e c)))
  659.      (setq y (n-term y))
  660.      (go a)
  661.     b     (unless (setq x (n-term x))
  662.         (return (pscheck (gvar-o xa) (list trunc) (cdr ans))))
  663.      (setq y yy a ans)
  664.     c     (when (or (null y) (e> (setq e (e+ (le x) (le y))) trunc))
  665.         (go b))
  666.      (setq c (pstimes (lc x) (lc y)))
  667.     d     (cond ((or (mono-term? a) (e> (le (n-term a)) e))
  668.         (add-term-&-pop a e c))
  669.            ((e> e (le (n-term a)))
  670.         (setq a (n-term a))
  671.         (go d))
  672.            (T (setq c (psplus c (lc (n-term a))))
  673.           (if (rczerop c) (rplacd a (n-term (n-term a)))
  674.              (change-lc (n-term a) c) (setq a (n-term a)))))
  675.      (setq y (n-term y))
  676.      (go c)))
  677.  
  678. (defun pscsubst (c v p)
  679.   (cond ((pscoefp p) p)
  680.     ((eq v (gvar p)) (pscsubst1 c p))
  681.     ((pointerp v (gvar p)) p)
  682.     (t (make-ps p (maplist
  683.                #'(lambda (q) (term (le q)
  684.                        (pscsubst c v (lc q))))
  685.                (terms p))))))
  686.  
  687. (defun pscsubst1 (v u)
  688.    (do ((a (rczero))
  689.     (ul (terms u) (n-term ul)))
  690.        ((null ul) a)
  691.       (setq a (psplus a (pstimes (lc ul) (psexpt v (le ul)))))))
  692.  
  693. (defun get-series (func trunc var e c)
  694.    (let ((pw (e// trunc e)))
  695.       (setq e (if (and (equal e (rcone)) (equal c (rcone)))
  696.           (getexp-fun func var pw)
  697.          (psmonsubst (getexp-fun func var pw) trunc e c)))
  698.       (if (and $float $keepfloat) (psfloat e) e)))
  699.  
  700. (defun psmonsubst (p trunc e c)
  701.   (if (psp p)
  702.       (psmonsubst1 p trunc e c
  703.            `(() . ,(terms p)) (ncons () ) (rcone) (rczero))
  704.     p))
  705.      
  706.  
  707. (defun psmonsubst1 (p trunc e c l ans cc el)
  708.    ;; We set $MAXTAYORDER to () here so that the calls to psexpt below
  709.    ;; won't do needless extra work, e.g. see rwg's complaint of 9/7/82.
  710.    (prog (a ee varh $maxtayorder)    
  711.      (setq a ans varh (gvar-o p))
  712.     a    (cond ((or (mono-term? l)
  713.             (e> (setq ee (e* e (le (n-term l)))) trunc))
  714.         (go end))
  715.            ((rczerop (setq cc
  716.                    (pstimes cc
  717.                     (psexpt c (e- (le (setq l (n-term l)))
  718.                               el))))))
  719.            ((mono-term? a)
  720.         (add-term a ee (pstimes cc (lc l)))))
  721.      (setq a (n-term a) el (le l))
  722.      (go a)
  723.     end  (return (pscheck varh (list trunc) (cdr ans)))))
  724.  
  725. (defun psexpon-gcd (terms)
  726.    (do ((gcd (le terms) (egcd (le l) gcd))
  727.     (l (n-term terms) (n-term l)))
  728.        ((null l) gcd)))
  729.  
  730. (defun psfind-s (p)
  731.    (if (psp p) (psfind-s (psterm (terms p) (rczero)))
  732.       (psfind-s1 p)))
  733.  
  734. (defun psfind-s1 (r)
  735.    (cond ((null (atom (cdr r))) (rczero))
  736.      ((atom (car r)) r)
  737.      (T (do ((p (pterm (cdar r) 0) (pterm (cdr p) 0)))
  738.         ((atom p) (cons p (cdr r)))))))
  739.  
  740. (defun psexpt (p n)
  741.     (cond ((rczerop n)            ;; p^0
  742.        (if (rczerop p)        ;; 0^0
  743.            (merror "~&Indeterminate form 0^0 generated inside PSEXPT~%")
  744.           (rcone)))            ;; Otherwise can let p^0 = 1
  745.       ((or (equal n (rcone)) (equal n (rcfone))) p)    ;; p^1 cases
  746.       ((pscoefp p) (rcexpt p n))
  747.       ((mono-term? (terms p))    ;; A monomial to a power
  748.        (let ((s (psfind-s n)) (n-s) (x) (l (terms p)))
  749.           ;; s is the numeric part of the exponent
  750.           (if (floatp (car s)) ;; Perhaps we souldn't
  751.           ;; rationalize if $keepfloat is true?
  752.           (setq s (MAXIMA-RATIONALIZE (*quo (car s) (cdr s)))))
  753.           (setq n-s (psdiff n s)    ;; the non-numeric part of exponent
  754.             x   (e* s (le l)))    ;; the degree of the lowest term
  755.           (setq x (if (and (null $maxtayorder) ;; if not getting all terms
  756.                    (e> x (t-o-var (gvar p))))
  757.               ;; and result is of high order
  758.               (rczero)    ;; then zero is enough
  759.              (pscheck (gvar-o p)    ;; otherwise
  760.                   (ncons (e+ (trunc-lvl p) ;; new trunc-level
  761.                          (e- x (le l)))) ;; kick exponent
  762.                   (ncons (term x (psexpt (lc l) n))))))
  763.           ;; x is now p^s
  764.           (if (or (rczerop n-s) (rczerop x))    ;; is that good enough?
  765.           x            ;; yes! The rest is bletcherous.
  766.          (pstimes x (psexpt (prep1 (m^ (get-inverse (gvar p))
  767.                            (rcdisrep n-s)))
  768.                     (ps-le p))))))
  769.       (T (prog (l lc le inc trunc s lt mr lim lcinv ans)
  770.            (setq lc (lc (setq l (terms p)))
  771.              le (le l) lt (lt l) trunc (trunc-lvl p)
  772.              inc (psexpon-gcd l) s (psfind-s n))
  773.            (when (floatp (car s))
  774.               (setq s (MAXIMA-RATIONALIZE (*quo (car s) (cdr s)))))
  775.            (setq ans (psexpt (setq lt (pscheck (gvar-o p) (list trunc)
  776.                                (list lt))) n)
  777.              lcinv (psexpt lc (rcmone))
  778.              mr (e+ inc (e* s le))
  779.              lim (if (and (infp trunc) (not (e> s (rczero))))
  780.                  (t-o-var (gvar p))
  781.                 ;; See the comment in PSEXPT1 below which tells
  782.                 ;; why we don't allow inf. trunc's here.
  783.                 (e+ (if (and (infp trunc) (not (rcintegerp s)))
  784.                     (if (infp (setq lim (t-o-var (gvar p))))
  785.                         (infin-ord-err)
  786.                        lim)
  787.                        trunc)
  788.                     (e* (e1- s) le)))
  789.              ans
  790.              (if (or (pscoefp ans) (null (eq (gvar p) (gvar ans))))
  791.                  (list 0 (term (rczero) ans))
  792.                 (cons 0 (terms ans))))
  793.            (and (null $maxtayorder)
  794.             (or (not (infp lim))
  795.                 (not (rcintegerp s))
  796.                 (e> (e* s (le (last l))) (t-o-var (gvar p))))
  797.             (setq lim (emin lim (t-o-var (gvar p)))))
  798.            ;;(and (infp lim) (n-term l) (e> (rczero) n)
  799.            ;;      (infin-ord-err))
  800.            (return (psexpt1 (gvar-o p)
  801.                     lim l n s inc 1 mr ans le lcinv))))))
  802.  
  803. (defun psexpt1 (varh trunc l n s inc m mr ans r linv)
  804.    ;; n is the power we are raising the series to
  805.    ;; inc is the exponent increment
  806.    ;; mr is the current exponent
  807.    ;; tr is the truncation level desired
  808.    ;; m is the term index
  809.    (declare (fixnum m))
  810.    ;; s ;Ignored <- not true, see below. Who wrote this?
  811.    (prog (a (k 0) ak cm-k c ma0 sum kr tr)
  812.      (declare (fixnum k))
  813.      ;; truly unfortunate that we need so many variables in this hack
  814.      (setq a (last ans) tr trunc)
  815.      ;; I don't see what's wrong with truncating exact series when
  816.      ;; raising them to fractional powers so we'll allow it for now.
  817.      ;; This is accomplished above in PSEXPT (see the comment). Thus,
  818.      ;; presumably, this check should never be needed anymore.
  819.      ;; Bugs catching this clause were sqrt(1-x)*taylor(f1,x,0,0)
  820.      ;; and sqrt(taylor(x+x^2,x,0,2)),taylor_truncate_polynomials=false.
  821.      (when (infp tr)
  822.         (if (rcintegerp s)
  823.         (setq tr (e* s (le (last l))))
  824.            (merror "Bad power series arg in PSEXPT")))
  825.      (when (infp tr) (setq tr (t-o-var (car varh))))
  826.      b (and (e> mr tr) (go end))
  827.        (setq kr inc ak l ma0 (pstimes (cons 1 m) linv)
  828.          k 1 sum (rczero))
  829.      a (if (or (> k m) (null (setq cm-k (psterm (cdr ans) (e- mr kr)))))
  830.            (go add-term))
  831.        (setq ak (or (pstrim-terms ak (e+ kr r)) (go add-term))
  832.          c (pstimes (psdiff (pstimes (cons k 1) n)
  833.                     (cons (f- m k) 1))
  834.                 (pstimes (if (e= (e+ kr r) (le ak))
  835.                      (lc ak)
  836.                      (rczero))
  837.                      cm-k)))
  838.        (setq sum (psplus sum c)
  839.          k (f1+ k) kr (e+ kr inc))
  840.        (go a)
  841.      add-term 
  842.       (and (null (rczerop sum))
  843.            (add-term-&-pop a mr (pstimes ma0 sum)))
  844.       (setq m (f1+ m) mr (e+ mr inc))
  845.      (go b)
  846.      end (return (pscheck varh (list trunc) (cdr ans)))))
  847.  
  848. (defun psderivative (p v)
  849.    (cond ((pscoefp p) (rcderiv p v))
  850.      ((eq v (gvar p))
  851.       (if (prog1 (rczerop (ps-le p))
  852.              (setq p (psderiv1 (gvar-o p)
  853.                 (trunc-lvl p) (cons 0 (terms p)) (list 0))))
  854.           (strip-zeroes p 'T) p))
  855.      (T (psderiv2 (gvar-o p)
  856.               (trunc-lvl p) v (cons 0 (terms p)) (list 0)))))
  857.  
  858. (defun psderiv1 (varh trunc l ans)
  859.        (do ((a (last ans)))
  860.        ((or (mono-term? l) (e> (le (n-term l)) trunc))
  861.         (pscheck varh (list (e1- trunc)) (cdr ans)))
  862.        (setq l (n-term l))
  863.        (when (not (rczerop (le l)))
  864.           (add-term-&-pop a (e1- (le l)) (pstimes (le l) (lc l))))))
  865.  
  866. (defun psderiv2 (varh trunc v l ans)
  867.        (do ((a (last ans) (n-term a)) (c))
  868.        ((or (mono-term? l) (e> (le (n-term l)) trunc))
  869.         (pscheck varh (list trunc) (cdr ans)))
  870.        (setq l (n-term l))
  871.        (or (rczerop (setq c (psderivative (lc l) v)))
  872.            (add-term a (le l) c))))
  873.  
  874. (defun psdp (p)
  875.   (let (temp temp2)
  876.    (cond ((pscoefp p) (rcderivx p))
  877.      ((or (rczerop (setq temp (getdiff (gvar-o p))))
  878.           (eq (car temp) 'MULTI))
  879.       (setq temp2 (psdp2 (gvar-o p) (trunc-lvl p)
  880.                  (cons 0 (terms p)) (list 0)))
  881.       (if (eq (car temp) 'MULTI)
  882.           (pstimes temp2
  883.                (make-ps (gvar-o p) (ncons (inf))
  884.                 (list (term (cdr temp) (rcone)))))
  885.           temp2))
  886.      (t (psdp1 (gvar-o p)
  887.            (trunc-lvl p) (cons 0 (terms p))
  888.            (list 0) temp)))))
  889.  
  890. (defun psdp1 (varh trunc l ans dx)
  891.        (do ((a (last ans)) (c (rczero)))
  892.        ((or (mono-term? l) (e> (le (n-term l)) trunc))
  893.         (psplus c (pscheck varh (list (e1- trunc)) (cdr ans))))
  894.        (setq l (n-term l))
  895.        (if (rczerop (le l)) (setq c (psdp (lc l)))
  896.            (add-term-&-pop
  897.         a (e1- (le l)) (pstimes (le l) (pstimes dx (lc l)))))))
  898.  
  899. (defun psdp2 (varh trunc l ans)
  900.        (do ((a (last ans)) (c))
  901.        ((or (mono-term? l) (e> (le (n-term l)) trunc))
  902.         (pscheck varh (list trunc) (cdr ans)))
  903.        (setq l (n-term l))
  904.        (when (null (rczerop (setq c (psdp (lc l)))))
  905.          (add-term-&-pop a (le l) c))))
  906.  
  907. ;;; Currently unused
  908. ;;;
  909. ;;; (defun psintegrate (p v)
  910. ;;;    (cond ((rczerop p) (rczero))
  911. ;;;      ((pscoefp p)
  912. ;;;       (pstimes p (taylor2 (get-inverse (car v)))))
  913. ;;;      ((eqgvar v (gvar-o p))
  914. ;;;       (psinteg1 (gvar-o p)
  915. ;;;             (trunc-lvl p) (cons 0 (terms p)) (list 0)))
  916. ;;;      (t (psinteg2 (gvar-o p)
  917. ;;;               (trunc-lvl p) v (cons 0 (terms p)) (list 0)))))
  918. ;;;
  919. ;;; (defun psinteg1 (varh trunc l ans)
  920. ;;;       (prog (a)
  921. ;;;         (setq a (last ans))
  922. ;;;    a    (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
  923. ;;;         (go end)
  924. ;;;         (add-term a (e1+ (le (setq l (n-term l))))
  925. ;;;               (pstimes (le l)
  926. ;;;                    (if (e= (le l) (rcmone))
  927. ;;;                    (prep1 (list '(%LOG)
  928. ;;;                             (get-inverse
  929. ;;;                              (car varh))))
  930. ;;;                    (lc l))))
  931. ;;;         (setq a (n-term a)))
  932. ;;;         (go a)
  933. ;;;        end  (return (pscheck varh (list (e1+ trunc)) (cdr ans)))))
  934.  
  935. ;;; (defun psinteg2 (varh trunc v l ans)
  936. ;;;        (prog (a)
  937. ;;;          (setq a (last ans))
  938. ;;;     a    (if (or (null (n-term l)) (e> (le (n-term l)) trunc))
  939. ;;;         (go end)
  940. ;;;         (add-term a (le l)
  941. ;;;               (psintegrate (lc (setq l (n-term l))) v))
  942. ;;;         (setq a (n-term a)))
  943. ;;;         (go a)
  944. ;;;    end  (return (pscheck varh (list trunc) (cdr ans)))))
  945.  
  946. (defun psexpt-log-ord (p)
  947.    (cond ((null $maxtayorder) (emin (trunc-lvl p) (t-o-var (gvar p))))
  948.      ((infp (trunc-lvl p)) (t-o-var (gvar p)))
  949.      (T (trunc-lvl p))))
  950.  
  951. (defun ps-infp (p)
  952.    (if (pscoefp p) ()
  953.       (get-) "..."))
  954.  
  955. (defun psexpt-fn (p)
  956.   (let (ans ord<0?)
  957.    (cond ((pscoefp p) (psexpt-fn2 (rcdisrep p)))
  958.      ((ps-lim-infp p) (psexpt-fn-sing p))
  959.      ((prog2 (setq ord<0? (e> (rczero) (ps-le p)))
  960.          (null (n-term (terms p))))
  961.       (setq ans (get-series '%EX (psexpt-log-ord p) (gvar-o p)
  962.                 (if ord<0? (e- (ps-le p)) (ps-le p))
  963.                 (ps-lc p)))
  964.       (if ord<0? (ps-invert-var ans) ans))
  965.      ((if ord<0?
  966.           (when (e= (rczero) (e (setq ans (ps-gt p))))
  967.          (pstimes (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p))
  968.                           (delq ans (terms p))))
  969.               (psexpt-fn2 (srdis (c ans)))))
  970.          (when (e= (rczero) (ps-le p))
  971.         (pstimes (psexpt-fn2 (srdis (lc (terms p))))
  972.              (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p))
  973.                          (n-term (terms p))))))) )
  974.      (T (prog (l le inc trunc lt ea0 ans)
  975.            (setq l (terms p))
  976.            (when ord<0?
  977.           ;(return (ps-invert-var (psexpt-fn (ps-invert-var p))))
  978.           (setq l (invert-terms l)))
  979.            (setq le (le l) lt (lt l) trunc (trunc-lvl p)
  980.              inc (psexpon-gcd l) ea0 (rcone))
  981.            (unless (e> (le l) (rczero))
  982.           (merror "Unreachable point in psexpt-fn"))
  983.            (setq ans
  984.              (if (or (pscoefp ea0) (null (eq (gvar p) (gvar ea0))))
  985.              (list 0 (term (rczero) ea0))
  986.             (cons 0 (terms ea0))))
  987.            (unless $maxtayorder
  988.           (setq trunc (emin trunc (t-o-var (gvar p)))))
  989.            (when (infp trunc) (setq trunc (t-o-var (gvar p))))
  990.            (setq ans (psexpt-fn1 (gvar-o p) trunc l inc 1 inc ans))
  991.            (return (if ord<0? (ps-invert-var ans) ans)))))))
  992.  
  993. (defun psexpt-fn-sing (p)
  994.    (let ((inf-var? (memq (gvar-lim (gvar p)) '($INF $MINF)))
  995.      (c*logs (c*logs (lt-poly p))) c strongest-term)
  996.       ;; Must pull out out logs here: exp(ci*log(ui)+x) -> ui^ci*exp(x)
  997.       ;; since its much harder for adjoin-tvar to do this transformation
  998.       ;; below after things have been disrepped.
  999.       (setq c (exp-c*logs c*logs) p (psdiff p (sum-c*logs c*logs)))
  1000.       (if (not (ps-lim-infp p))
  1001.       ;; Here we just subtracted the only infinite term, e.g.
  1002.       ;; p = 1/2*log(x)+1/log(x)+...
  1003.       (pstimes c (psexpt-fn p))
  1004.      (setq strongest-term (if inf-var? (ps-gt p) (ps-lt p)))
  1005.      ;; If the strongest term has degree 0 in the mainvar then the singular
  1006.      ;; terms occur in some other weaker var. There may be terms in this
  1007.      ;; coef which arent singular (e.g. 1 in (1/x+1+...)+exp(-1/x)+...) so
  1008.      ;; we must recursively psexpt-fn this term to get only what we need.
  1009.      (if (rczerop (e strongest-term))
  1010.          (setq c (pstimes c (psexpt-fn (c strongest-term))))
  1011.         (dolist (exp (expand-and-disrep strongest-term p))
  1012.            (setq c (pstimes c (adjoin-tvar (m^ '$%E exp))))))
  1013.      (pstimes c (psexpt-fn (pscheck (gvar-o p) (list (trunc-lvl p))
  1014.                     (if inf-var?
  1015.                         (delq strongest-term (terms p))
  1016.                        (n-term (terms p)))))))))
  1017.  
  1018. (defun gvar-logp (gvar)
  1019.    (let ((var (get-inverse gvar)))
  1020.       (and (consp var) (eq (caar var) 'MEXPT) (equal (caddr var) -1)
  1021.        (consp (setq var (cadr var))) (eq (caar var) '%LOG)
  1022.        var)))
  1023.  
  1024. (defun c*logs (p)
  1025.    (if (pscoefp p) ()
  1026.       (let ((log (gvar-logp (gvar p))) c)
  1027.      (if (not log) ()
  1028.         (setq c (psconst (psterm (terms p) (rcmone))))
  1029.         ;; We don't want e.g. exp(x^a*log(x)) -> x^x^a
  1030.         (if (not (mfree (rcdisrep c) tvars)) ()
  1031.            (cons (cons c (cons log p))
  1032.              (c*logs (psterm (terms p) (rczero)))))))))
  1033.  
  1034. (defun psconst (p)
  1035.    (if (pscoefp p) p (psconst (psterm (terms p) (rczero)))))
  1036.  
  1037. (defun exp-c*logs (c*logs)
  1038.    (if (null c*logs) (rcone)
  1039.       (pstimes (taylor2 `((MEXPT) ,(cadr (cadr (car c*logs)))
  1040.                   ,(rcdisrep (caar c*logs))))
  1041.            (exp-c*logs (cdr c*logs)))))
  1042.  
  1043. (defun sum-c*logs (c*logs)
  1044.    (if (null c*logs) (rczero)
  1045.       (let ((ps (cddr (car c*logs))))
  1046.      (psplus (make-ps ps (ncons (term (rcmone) (caar c*logs))))
  1047.          (sum-c*logs (cdr c*logs))))))
  1048. #|
  1049. (defun remove-c*logs (p)
  1050.    (if (pscoefp p) (rcone)
  1051.       (let ((log (gvar-logp (gvar p)))
  1052.         (c1 (psterm (terms p) (rcone))) const)
  1053.      (if (pscoefp c1)
  1054.          (if (not (mfree (rcdisrep c0) tvars)) (rczero)
  1055.         (setq p (remterm p (rcone)))
  1056.         c1)
  1057.         (pstimes (if log (taylor2 `((MEXPT) ,(cadr log)
  1058.                         ,(rempsconst
  1059.                           (psterm (terms p) (rcone)))))
  1060.             (rcone))
  1061.              (remove-c*logs (psterm (terms p) (rczero))))))))
  1062.  
  1063. (defun rempsconst (ps)
  1064.    (let ((c0 (psterm (terms ps) (rczero))))
  1065.       (if (psp c0) (rempsconst c0)
  1066.      (if (not (mfree (rcdisrep c0) tvars)) (rczero)
  1067.         (remterm (terms ps) (rczero))
  1068.         c0))))
  1069.  
  1070. (defun remterm (terms e)
  1071.    (do ((terms* terms (n-term terms*)))
  1072.        ((null terms) (merror "remterm ran out of terms"))
  1073.       (when (equal (le terms*) e)
  1074.      (return (setq terms (delq (lt terms*) terms))))))
  1075. |#
  1076. ;; Calculatest the limit of a series at the expansion point. Returns one of
  1077. ;; {$ZEROA, $ZEROB, $POS, $NEG, $INF, $MINF}.
  1078.  
  1079. (defvar tvar-limits ()
  1080.    "A list of the form ((gvar . limit(gvar)) ...)")
  1081. #|
  1082. (defun ps-lim (ps)
  1083.    (if (pscoefp ps) (coef-sign (srdis ps))
  1084.       ;; Assume taylor vars at 0+ for now. Should handle the cases when
  1085.       ;; the expansion point is INF, MINF,etc.
  1086.       (let* ((lim (gvar-lim (gvar ps)))
  1087.          (strongest-term
  1088.           (if (memq lim '($INF $MINF)) (ps-gt ps) (ps-lt ps))))
  1089.      (if (ezerop (e strongest-term)) (ps-lim (c strongest-term))
  1090.         (setq lim (lim-power lim (e strongest-term)))
  1091.         ;; It is assumed that gvars of coeff's are weaker than the gvar
  1092.         ;; of a ps so they only affect the sign of the overall limit.
  1093.         (if (lim-plusp (ps-lim (c strongest-term))) lim
  1094.            (lim-minus lim)))))) |#
  1095.  
  1096. (defun ps-lim-infp (ps)
  1097.    (if (pscoefp ps) ()
  1098.       ;; Assume taylor vars at 0+ for now. Should handle the cases when
  1099.       ;; the expansion point is INF, MINF,etc.
  1100.       (let* ((lim (gvar-lim (gvar ps)))
  1101.          (strongest-term
  1102.           (if (memq lim '($INF $MINF)) (ps-gt ps) (ps-lt ps))))
  1103.      (if (ezerop (e strongest-term)) (ps-lim-infp (c strongest-term))
  1104.         (setq lim (lim-power lim (e strongest-term)))
  1105.         (and (lim-infp lim) (not (eq lim '$INFINITY)))))))
  1106.  
  1107. (defun lim-zerop (lim) (memq lim '($ZEROA $ZEROB $ZEROIM)))
  1108. (defun lim-plusp (lim) (memq lim '($ZEROA $POS $INF $FINITE)))
  1109. (defun lim-finitep (lim) (memq lim '($POS $NEG $IM $FINITE))) 
  1110. (defun lim-infp (lim) (memq lim '($INF $MINF $INFINITY)))
  1111. (defun lim-imagp (lim) (memq lim '($IM $INFINITY)))
  1112. (defun lim-minus (lim)
  1113.    (cdr (assq lim '(($zeroa . $zerob) ($zerob . $zeroa) ($pos . $neg)
  1114.             ($neg . $pos) ($inf . $minf) ($minf . $inf)
  1115.             ($im . $im) ($infinity . $infinity) ($finite . $finite)))))
  1116. (defun lim-abs (lim)
  1117.    (or (cdr (assq lim '(($zerob . $zeroa) ($neg . $pos) ($minf . $inf))))
  1118.        lim))
  1119.  
  1120. (defun lim-times (lim1 lim2)
  1121.   (let (lim) 
  1122.    (cond ((and (lim-infp lim1) (lim-infp lim2)) (setq lim '$INF))
  1123.      ((and (lim-zerop lim1) (lim-zerop lim2)) (setq lim '$POS))
  1124.      ((or (when (lim-finitep lim2) (exch lim1 lim2) 'T)
  1125.           (lim-finitep lim1))
  1126.       (when (and (eq lim1 '$finite) (lim-infp lim1))
  1127.          (break "Undefined finite*inf in lim-times")) 
  1128.       (setq lim (lim-abs lim2)))
  1129.      (T (break "Undefined limit product ~A * ~A in lim-times lim1 lim2")))
  1130.    (if (or (lim-imagp lim1) (lim-imagp lim2))
  1131.        (if (lim-infp lim) '$INFINITY '$IM)
  1132.       (if (and (lim-plusp lim1) (lim-plusp lim2)) lim (lim-minus lim)))))
  1133.  
  1134. (defun lim-power (lim power)
  1135.    (cond ((ezerop power) '$POS)
  1136.      ((e> (rczero) power) (lim-recip (lim-power lim (e- power))))
  1137.      ((not (oddp (car power)))
  1138.       (if (lim-plusp lim) lim (lim-minus lim)))
  1139.      (T lim)))
  1140.  
  1141. (defun lim-recip (lim)
  1142.    (or (cdr (assq lim '(($ZEROA . $INF) ($ZEROB . $MINF)
  1143.                 ($INF . $ZEROA) ($MINF . $ZEROB))))
  1144.        (if (eq lim '$finite) (break "inverting $FINITE?")
  1145.       lim)))
  1146.  
  1147. (defun lim-exp (lim)
  1148.    (case lim
  1149.       (($ZEROA $ZEROB $POS $NEG $MINF) '$ZEROA)
  1150.       (($INF $FINITE) lim)
  1151.       ($INFINITY '$INFINITY) ; actually only if Re lim = inf
  1152.       (T (break "Unhandled limit in lim-exp"))))
  1153.  
  1154. (defun lim-log (lim)
  1155.    (case lim
  1156.       ($ZEROA '$MINF)
  1157.       ($INF '$INF)
  1158.       (T (break "Unhandled limit in lim-log"))))
  1159.       
  1160. (defun expand-and-disrep (term p)
  1161.    (let ((x^n (list '(MEXPT) (get-inverse (gvar p)) (edisrep (e term))))
  1162.      (a (c term)))
  1163.       (if (pscoefp a) (ncons (m* (srdis a) x^n))
  1164.      (mapcar #'(lambda (subterm)
  1165.               (m* (cons '(MTIMES) (expand-and-disrep subterm a)) x^n))
  1166.          (terms a)))))
  1167.  
  1168. (defun adjoin-sing-datum (d)
  1169.    (let ((r (prep1 (datum-var d))) g* (g (gensym)) (kernel (datum-var d))
  1170.          (no (f1+ (cdr (int-var (car (last tlist)))))))
  1171.       (unless (and (equal (car r) 1) (equal (cddr r) '(1 1)))
  1172.      (break "bad singular datum"))
  1173.       (putprop g kernel 'disrep)
  1174.       (rplacd (cdddr d) (cons g no))
  1175.       (adjoin-datum d)
  1176.       (push (cons (setq g* (cadr r)) kernel) key-vars)
  1177.       (push (cons g kernel) key-vars)
  1178.       (push (car key-vars) ivars)
  1179.       ;(push (cons kernel (cons (pget g) 1)) genpairs)
  1180.       (push (cons g (exp-pt d)) tvar-limits)))   
  1181.  
  1182. (defun adjoin-tvar (exp) (rat->ps (prep1 exp)))
  1183.  
  1184. (defun rat->ps (rat)
  1185.    (pstimes (poly->ps (car rat))
  1186.         (psexpt (poly->ps (cdr rat)) (rcmone))))
  1187.  
  1188. (defun poly->ps (poly)
  1189.    (if (or (pcoefp poly) (mfree (pdis poly) tvars)) (prep1 poly)
  1190.       (let ((g (p-var poly)) datum (pow (rcone)))
  1191.      (if (setq datum (key-var-pow g)) (desetq (g . pow) datum)
  1192.         (desetq (g . pow) (adjoin-pvar g)))
  1193.     (if (and (not (atom g)) (psp g)) g
  1194.      (setq datum (gvar-data g))
  1195.      (do ((po-terms (p-terms poly) (p-red po-terms))
  1196.           (ps-terms () (push (term (e* pow (prep1 (pt-le po-terms)))
  1197.                        (poly->ps (pt-lc po-terms)))
  1198.                  ps-terms)))
  1199.          ((null po-terms)
  1200.           ;; This must be exact so that when we invert in rat-ps above
  1201.           ;; we dont lose any terms. E.g. try
  1202.           ;; taylor(log(1+exp(-1/x)),x,0,5). When taylor2'ing exp(-1/x),
  1203.           ;; if you used current trunc here this would return exp(1/x)...5
  1204.           ;; which would then be trunc'd to degree 3 by psexpt.
  1205.           (make-ps (int-var datum)
  1206.                (ncons (current-trunc datum))
  1207.                (if (eq g (data-gvar datum)) ps-terms
  1208.               (invert-terms ps-terms)))))))))
  1209.  
  1210. (defun key-var-pow (g)
  1211.    (let ((var (get-key-var g)) datum)
  1212.       (when var
  1213.      (setq datum (get-datum var))
  1214.      (if (eq g (setq g (data-gvar datum))) (cons g (rcone))
  1215.         (cons g (rcmone))))))
  1216.  
  1217. (defun adjoin-pvar (g)
  1218.    (let ((kernel (get g 'disrep)) g* lim datum ans
  1219.      (no (f1+ (cdr (int-var (car (last tlist)))))) (pow (rcone)) expt)
  1220.       (when (assol kernel tlist) (break 'bad1))
  1221.      (if (and (eq (caar kernel) 'MEXPT) (eq (cadr kernel) '$%E)
  1222.           (not (atom (setq expt (caddr kernel))))
  1223.           (eq (caar expt) 'MTIMES) (not (mfree expt (ncons '$%I))))
  1224.      (if (not (and (eq (cadr expt) '$%I) (mfree (cddr expt) (ncons '$%I))))
  1225.          (break "bad exponential in adjoin-pvar") (break "now")
  1226.         (setq expt (m// expt '$%I))
  1227.         (cons (psplus (adjoin-tvar `((%COS) ,expt))
  1228.               (pstimes (prep1 '$%I) (adjoin-tvar `((%SIN) ,expt))))
  1229.           pow))
  1230.       (when (eq (caar kernel) 'MEXPT)
  1231.      (when (and (not (atom (setq expt (caddr kernel))))
  1232.             (eq (caar expt) 'MTIMES)
  1233.             ($ratnump (cadr expt)))
  1234.         (setq pow (cadr expt) kernel (m^ kernel (m// pow))
  1235.           g (prep1 kernel) pow (prep1 pow))
  1236.         (unless (and (equal (cdr g) 1) (equal (cdar g) '(1 1)))
  1237.            (break "Illegal kernel in ADJOIN-PVAR"))
  1238.         (setq g (caar g) kernel (get g 'disrep))))
  1239.       (if (setq ans (key-var-pow g)) (cons (car ans) (e* pow (cdr ans)))
  1240.      (when (lim-infp (or lim (setq lim (tvar-lim kernel))))
  1241.         (setq g* g g (gensym) kernel (m// kernel)
  1242.           lim (lim-recip lim) pow (e* (rcmone) pow))
  1243.         (putprop g kernel 'disrep)
  1244.         ;(push g genvar) (push kernel varlist)
  1245.         (push (cons g* kernel) key-vars))
  1246.      (when (assol kernel tlist) (break 'bad2))
  1247.      (setq datum (list* kernel
  1248.                 ;(mapcar #'(lambda (e) (emax e (rczero)))
  1249.                 ;        (trunc-stack (car tlist)))
  1250.                 (copy-top-level (trunc-stack (car tlist)))
  1251.                 lim () g no))
  1252.      ;(setq tlist (nconc tlist (ncons datum)))
  1253.      (adjoin-datum datum)
  1254.      (push (cons g kernel) key-vars)
  1255.      (push (car key-vars) ivars)
  1256.      ;(push (cons kernel (cons (pget g) 1)) genpairs)
  1257.      (push (cons g lim) tvar-limits)
  1258.      (cons g pow)))))
  1259.  
  1260. (defun adjoin-datum (datum)
  1261.    (do ((tlist* tlist (cdr tlist*))
  1262.     (tlist** () tlist*))
  1263.        ((null tlist*) (setq tlist (nconc tlist (ncons datum))))
  1264.       (when (stronger-var? (datum-var (car tlist*)) (datum-var datum))
  1265.      (return (if (null tlist**)
  1266.              (progn (push datum tlist) (renumber-tlist tlist))
  1267.             (rplacd tlist** (cons datum tlist*))
  1268.             (renumber-tlist (cdr tlist**)))))))
  1269.  
  1270. ;; Maybe this should just permute the numbering in case it isn't sequential?
  1271.  
  1272. (defun renumber-tlist (tlist)
  1273.    (rplacd (data-gvar-o (car tlist)) (cdr (data-gvar-o (cadr tlist))))
  1274.    (do ((tlist* (cdr tlist) (cdr tlist*)))
  1275.        ((null tlist*))
  1276.       (rplacd (data-gvar-o (car tlist*))
  1277.           (f1+ (cdr (data-gvar-o (car tlist*)))))))
  1278.  
  1279. (defun tvar? (var) (or (atom var) (memq 'array (cdar var))))
  1280.  
  1281. ;; Needs to be extended to handle singular tvars in > 1 var's.
  1282.  
  1283. (defun stronger-var? (v1 v2)
  1284.    (let ((e1 (rcone)) (e2 (rcone)) reverse? ans)
  1285.       (when (alike1 v1 v2) (break "stronger-var? called on equal vars"))
  1286.       (when (and (mexptp v1) ($ratnump (caddr v1)))
  1287.      (setq e1 (prep1 (caddr v1)) v1 (cadr v1)))
  1288.       (when (and (mexptp v2) ($ratnump (caddr v2)))
  1289.      (setq e2 (prep1 (caddr v2)) v2 (cadr v2)))
  1290.     (if (alike1 v1 v2)
  1291.     (if (equal e1 e2) (break "= vars in stronger-var?")
  1292.        (e> e1 e2))
  1293.      (when (eq (tvar-lim v2) '$FINITE)
  1294.     (exch v1 v2) (exch e1 e2) (setq reverse? (not reverse?)))
  1295.     (if (eq (tvar-lim v1) '$FINITE)
  1296.     (if (eq (tvar-lim v2) '$FINITE) (great v1 v2) reverse?)
  1297.      (when (mtimesp v2)
  1298.      (exch v1 v2) (exch e1 e2) (setq reverse? (not reverse?)))
  1299.      (setq ans
  1300.       (if (mtimesp v1)
  1301.       (stronger-vars? (order-vars-by-strength (cdr v1))
  1302.          (order-vars-by-strength (if (mtimesp v2) (cdr v2)
  1303.                     (ncons (m^ v2 (edisrep e2))))))
  1304.      (when (tvar? v2)
  1305.         (exch v1 v2) (exch e1 e2) (setq reverse? (not reverse?)))
  1306.        (if (tvar? v1)
  1307.       (cond ((tvar? v2)
  1308.          (let ((n1 (cdr (data-gvar-o (get-datum v1 T))))
  1309.                (n2 (cdr (data-gvar-o (get-datum v2 T)))))
  1310.             (> n1 n2)))
  1311.         ((mfree v2 (ncons v1))
  1312.          (break "Unhandled multivar datum comparison"))
  1313.         ((eq (caar v2) '%LOG) 'T)
  1314.         ((and (eq (caar v2) 'MEXPT) (eq (cadr v2) '$%E))
  1315.          (stronger-var? `((%LOG) ,v1) (caddr v2)))
  1316.         (T (break "Unhandled var in stronger-var?")))
  1317.      (when (eq (caar v2) '%LOG)
  1318.         (exch v1 v2) (exch e1 e2) (setq reverse? (not reverse?)))
  1319.      (if (eq (caar v1) '%LOG)
  1320.          (cond ((eq (caar v2) '%LOG)
  1321.             (stronger-var? (log-abs-tvar v1) (log-abs-tvar v2)))
  1322.            ((and (eq (caar v2) 'MEXPT) (eq (cadr v2) '$%E))
  1323.             (stronger-var? `((%LOG) ,v1) (caddr v2)))
  1324.            (T (break "Unhandled var in stronger-var?")))
  1325.         (if (and (eq (caar v1) 'MEXPT) (eq (cadr v1) '$%E)
  1326.              (eq (caar v2) 'MEXPT) (eq (cadr v2) '$%E))
  1327.         (stronger-var? (caddr v1) (caddr v2))
  1328.            (break "Unhandled var in stronger-var?"))))))
  1329.      (if reverse? (not ans) ans)))))
  1330.  
  1331. (defun neg-monom? (exp)
  1332.    (and (mtimesp exp) (equal (cadr exp) -1) (null (cdddr exp))
  1333.     (caddr exp)))
  1334.  
  1335. (defun log-abs-tvar (var)
  1336.    (cond ((or (tvar? (cadr var)) (eq (caar (cadr var)) '%LOG)) (cadr var))
  1337.      ((neg-monom? (cadr var)) )
  1338.      (T (break "Illegal logarithmic tvar"))))
  1339.  
  1340. (defun order-vars-by-strength (vars)
  1341.    (do ((vars* vars (cdr vars*)) (ordvars () ))
  1342.        ((null vars*) ordvars)
  1343.       (unless (mfree (car vars*) tvars) ; ignore constants
  1344.      (do ((ordvars* ordvars (cdr ordvars*)))
  1345.          ((null ordvars*)
  1346.           (if (null ordvars) (setq ordvars (ncons (car vars*)))
  1347.          (rplacd (last ordvars) (ncons (car vars*)))))
  1348.         (when (stronger-var? (car vars*) (car ordvars*))
  1349.            (rplacd ordvars* (cons (car ordvars*) (cdr ordvars*)))
  1350.            (rplaca ordvars* (car vars*))
  1351.            (return () ))))))
  1352.  
  1353. (defun stronger-vars? (vars1 vars2)
  1354.    (do ((vars1* vars1 (cdr vars1*))
  1355.     (vars2* vars2 (cdr vars2*)))
  1356.        (())
  1357.       (cond ((null vars1*)
  1358.          (if (null vars2*) (break "two equal vars generated")
  1359.         (let ((lim (tvar-lim (car vars2*))))
  1360.            (return
  1361.             (cond ((lim-infp lim) ())
  1362.               ((lim-zerop lim) 'T)
  1363.               (T (break "var with non-zero finite lim?")))))))
  1364.         ((null vars2*)
  1365.          (if (null vars1*) (break "two equal vars generated")
  1366.         (let ((lim (tvar-lim (car vars1*))))
  1367.            (return
  1368.             (cond ((lim-infp lim) 'T)
  1369.               ((lim-zerop lim) ())
  1370.               (T (break "var with non-zero finite lim?")))))))
  1371.         ((alike1 (car vars1*) (car vars2*)) )
  1372.         ((return (stronger-var? (car vars1*) (car vars2*)))))))
  1373.  
  1374. (defun stronger-datum? (d1 d2)
  1375.    (setq d1 (datum-var d1) d2 (datum-var d2))
  1376.    (do ((end-flag) (answer))
  1377.        (end-flag (memq answer '($yes $y)))
  1378.       (setq answer (retrieve `((mtext) |Is  | ,d1 | stronger than | ,d2 |?|)
  1379.                  nil))
  1380.       (if (memq answer '($yes $y $no $n)) (setq end-flag 'T)
  1381.      (mtell "~%Acceptable answers are Yes, Y, No, N~%"))))
  1382.  
  1383. (defun datum-lim (datum)
  1384.    (if (not (tvar? (datum-var datum)))
  1385.        ;(cdr (assq (data-gvar datum) tvar-limits))
  1386.        (exp-pt datum)
  1387.       (let ((pt (exp-pt datum)))
  1388.      (if (memq pt '($INF $MINF)) pt '$ZEROA))))
  1389.  
  1390. (defun tvar-lim (kernel)
  1391.   (if (mfree kernel tvars) (coef-sign kernel)
  1392.     (let ((datum (get-datum kernel T)) lim)
  1393.       (or (and datum (datum-lim datum))
  1394.       (and (setq datum (get-datum (m// kernel) T))
  1395.            (setq lim (datum-lim datum))
  1396.            (lim-recip lim))
  1397.       (progn
  1398.        (setq lim
  1399.          (cond ((eq (caar kernel) 'MEXPT)
  1400.             (cond ((and (setq datum (get-datum (cadr kernel) T))
  1401.                     ($ratnump (caddr kernel)))
  1402.                    (lim-power (datum-lim datum)
  1403.                       (prep1 (caddr kernel))))
  1404.                   (($ratnump (caddr kernel))
  1405.                    (lim-power (tvar-lim (cadr kernel))
  1406.                       (prep1 (caddr kernel))))
  1407.                   ((eq (cadr kernel) '$%E)
  1408.                    (lim-exp (tvar-lim (caddr kernel))))
  1409.                   (T (break "Unhandled case in tvar-lim"))))
  1410.                ((eq (caar kernel) 'MTIMES)
  1411.             (do ((ans (tvar-lim (cadr kernel))
  1412.                   (lim-times ans (tvar-lim (car facs))))
  1413.                  (facs (cddr kernel) (cdr facs)))
  1414.                 ((null facs) ans)))
  1415.                ((eq (caar kernel) '%LOG)
  1416.             ;; Assume all log's are of the form log(x+a),
  1417.             ;; log(-log(x+a)),... First type go to minf; all
  1418.             ;; others to inf.
  1419.             ;(lim-log (tvar-lim (cadr kernel)))
  1420.             (if (tvar? (cadr kernel))
  1421.                 (let ((pt (exp-pt (get-datum (cadr kernel)))))
  1422.                    (if (memq pt '($INF $MINF)) (lim-log pt)
  1423.                   '$MINF))
  1424.                (cond ((eq (caar (cadr kernel)) 'MPLUS) '$MINF)
  1425.                  ((eq (caar (cadr kernel)) '%LOG) '$INF)
  1426.                  ((neg-monom? (cadr kernel)) '$INF)
  1427.                  (T (break "Illegal log kernel")))))
  1428.                ((memq (caar kernel) '(%SIN %COS))
  1429.             (unless (lim-infp (tvar-lim (cadr kernel)))
  1430.                (break "Invalid trig kernel in tvar-lim"))
  1431.             '$FINITE)
  1432.                (T (break "Unhandled kernel in tvar-lim"))))
  1433.        (when datum (push (cons (data-gvar datum) lim) tvar-limits))
  1434.        lim)))))
  1435.  
  1436. (defun coef-sign (coef)
  1437.    (if (not ($freeof '$%I coef)) '$IM ($asksign coef)))
  1438.  
  1439. (defun gvar-lim (gvar)
  1440.    (or (cdr (assq gvar tvar-limits))
  1441.        (if (memq (gvar->var gvar) tvars) '$ZEROA ; user tvars assumed 0+ now
  1442.       (break "Invalid gvar"))))
  1443.  
  1444. (defun psexpt-fn1 (varh trunc l inc m mr ans)
  1445.        (declare (fixnum m ))
  1446.        (prog (a (k 0) ak cm-k c sum kr lim)
  1447.          (declare (fixnum k ))
  1448.            ;; truly unfortunate that we need so many variables in this hack
  1449.        (setq a (last ans))
  1450.      b (and (e> mr trunc) (go end))
  1451.        (setq kr inc ak l k 1 sum (rczero) lim m)
  1452.      a (cond ((or (> k lim)
  1453.               (null (setq cm-k (psterm (cdr ans) (e- mr kr)))))
  1454.           (go add-term)))
  1455.        (setq ak (or (pstrim-terms ak kr)
  1456.             (go add-term))
  1457.          c (pstimes (ereduce k m)
  1458.                 (pstimes (psterm1 ak kr) cm-k))
  1459.          sum (psplus sum c))
  1460.        (setq k (f1+ k) kr (e+ kr inc))
  1461.        (go a)
  1462.      add-term 
  1463.        (unless (rczerop sum) (add-term-&-pop a mr sum))
  1464.        (setq m (f1+ m) mr (e+ mr inc))
  1465.        (go b)
  1466.      end
  1467.        (return (pscheck varh (list trunc) (cdr ans)))))
  1468.  
  1469. ;;; PSEXPT-FN2 and RED-MONO-LOG are needed to reduce exponentials of logs.
  1470.  
  1471. (defun psexpt-fn2 (p)
  1472.    (cond ((atom p) (if (get-datum p) (psexpt-fn (taylor2 p))
  1473.               (prep1 `((MEXPT) $%E ,p))))
  1474.      ((eq (caar p) '%LOG)
  1475.       (if (get-datum (cadr p)) (taylor2 (cadr p)) (prep1 (cadr p))))
  1476.      ((or (eq (caar p) 'MPLUS) (eq (caar p) 'MTIMES))
  1477.       (let ((e ($ratexpand p)) temp)
  1478.          (ifn (and (consp e) (memq (caar e) '(MPLUS MTIMES)))
  1479.           (psexpt-fn2 e) ;; handles things like '((MPLUS) $Y) --> $Y
  1480.         (if (eq (caar e) 'MPLUS)
  1481.             (do ((sumnds (cdr e) (cdr sumnds))
  1482.              (log-facs) (l))
  1483.             ((null sumnds)
  1484.              (ifn log-facs (tsexpt '$%E p)
  1485.                 (tstimes (cons (m^t '$%E (m+l l)) log-facs))))
  1486.                (if (setq temp (red-mono-log (car sumnds)))
  1487.                (push temp log-facs)
  1488.               (push (car sumnds) l)))
  1489.            (setq temp (red-mono-log e)) 
  1490.            (if temp (taylor2 temp) (prep1 (power '$%E p)))))))
  1491.      (T (prep1 (power '$%E p)))))
  1492.  
  1493. (defun red-mono-log (e)
  1494.    (cond ((atom e) ())
  1495.      ((eq (caar e) '%LOG) (cadr e))
  1496.      ((mtimesp e)
  1497.       (do ((facs (cdr e) (cdr facs)) (log-term))
  1498.           ((null facs)
  1499.            (when log-term
  1500.              (m^t (cadr log-term) (m*l (remq log-term (cdr e))))))
  1501.           (if (and (null (atom (car facs))) (eq (caaar facs) '%LOG))
  1502.           (if log-term (return ()) (setq log-term (car facs)))
  1503.           (ifn (mfree (car facs) tvars) (return () )))))
  1504.      (T () )))
  1505.  
  1506. (defun pslog (p)
  1507.    (if (pscoefp p) (pslog2 (rcdisrep p))
  1508.        (let ((terms (terms p)))
  1509.       (cond ((mono-term? terms) ; log(c x^n) = log(c) + n log(x)
  1510.          ;; do this always for now
  1511.          (if 'T ;$TAYLOR_LOGEXPAND
  1512.              ;(psplus (pslog (lc terms))
  1513.              ;         (pstimes (le terms) (pslog-of-gvar (gvar p))))
  1514.              (pslog-monom p)
  1515.              (prep1 `((%LOG) ,(term-disrep (lt terms) p)))))
  1516.         ;; expand log(1+ax^n) directly by series substitution
  1517.         ((not (or (n-term (setq terms (terms (psplus p (rcmone)))))
  1518.               ;(e> (rczero) (le terms))
  1519.               (ps-lim-infp p)))
  1520.          (setq p (get-series '%LOG (psexpt-log-ord p) (gvar-o p)
  1521.                 (if (e> (rczero) (le terms)) (e- (le terms))
  1522.                    (le terms))
  1523.                 (lc terms)))
  1524.          (if (e> (rczero) (le terms)) (ps-invert-var p) p))
  1525.         (T (prog (l inc trunc lt ans lterm $maxtayorder gvar-lim gt)
  1526.             ;; log(lt+y) = log(lt) + log(1 + y/lt) = lterm + p
  1527.             (setq trunc (trunc-lvl p))
  1528.             (if (not (memq (setq gvar-lim (gvar-lim (gvar p)))
  1529.                    '($ZEROA $ZEROB $INF $MINF)))
  1530.             (break "bad gvar lim")
  1531.                (if (memq gvar-lim '($INF $MINF))
  1532.                (setq lt (ps-gt p) gt lt)
  1533.               (setq lt (ps-lt p) gt () )))
  1534.             (setq lterm (pslog
  1535.                  (setq lt (pscheck (gvar-o p)
  1536.                            (ncons trunc)
  1537.                            (ncons lt))))
  1538.               p (pstimes p (let (($MAXTAYORDER 'T))
  1539.                       (psexpt lt (rcmone)))))
  1540.             (when (and (memq gvar-lim '($INF $MINF))
  1541.                    (e> (le terms) (rczero)))
  1542.                (return (psplus lterm (pslog p))))
  1543.             (when (pscoefp p)
  1544.                (unless (equal p (rcone))
  1545.               (merror "Internal TAYLOR MAXIMA-ERROR in PSLOG"))
  1546.                (return lterm))
  1547.             (setq l (terms p) inc (psexpon-gcd l))
  1548.             (if gt (setq l (zl-DELETE (last l) l))
  1549.                (setq l (n-term l)))
  1550.             (setq ans (ncons 0))
  1551.             (unless $maxtayorder
  1552.                (setq trunc (emin trunc (t-o-var (gvar p)))))
  1553.             ;; When we've divided by the greatest term, all terms
  1554.             ;; have non-positive exponents and we must perform the
  1555.             ;; transformation x -> 1/x befor calling pslog1 and then
  1556.             ;; perform the inverse afterwards.
  1557.             (when gt (setq l (invert-terms l)))
  1558.             (when (e> (rczero) inc) (setq inc (e- inc)))
  1559.             (setq ans (psplus lterm
  1560.                  (pslog1 (gvar-o p) trunc l inc 1 inc ans)))
  1561.             (return
  1562.              (if (and gt (psp ans) (eq (gvar ans) (gvar p)))
  1563.              (ps-invert-var ans)
  1564.             ans))))))))
  1565.  
  1566. (defun invert-terms (terms)
  1567.    (nreverse (mapc #'(lambda (x) (rplaca x (e- (e x)))) terms)))
  1568.  
  1569. (defun ps-invert-var (ps)
  1570.    (when (psp ps) (rplacd (cddr ps) (invert-terms (terms ps))))
  1571.    ps)
  1572.  
  1573. (defun ps-gt (ps)
  1574.    (if (pscoefp ps) (term (rczero) ps)
  1575.       (lt (last (terms ps)))))
  1576.  
  1577. (defun pslog1  (varh trunc l inc m mr ans)
  1578.        (declare (fixnum m ))
  1579.        (prog (a (k 0) ak cm-k c sum kr m-kr)
  1580.          (declare (fixnum k ))
  1581.        ;; truly unfortunate that we need so many variables in this hack
  1582.        ;; 
  1583.        (setq a (last ans))
  1584.      b (and (e> mr trunc) (go end))
  1585.        (setq kr inc ak l k 1 sum (rczero))
  1586.      a (cond ((or (= k m)
  1587.               (null (setq cm-k (psterm (cdr ans)
  1588.                            (setq m-kr (e- mr kr))))))
  1589.           (go add-term)))
  1590.        (setq ak (or (pstrim-terms ak kr)
  1591.             (go add-term))
  1592.          c (pstimes m-kr (pstimes (psterm1 ak kr) cm-k))
  1593.          sum (psplus sum c)
  1594.          k (f1+ k) kr (e+ kr inc))
  1595.        (go a)
  1596.      add-term 
  1597.        (cond ((setq c (pstrim-terms ak mr))
  1598.           (setq c (psterm1 c mr)))
  1599.          ((setq c (rczero))))        
  1600.        (setq sum (psdiff c (pstimes sum (e// mr))))
  1601.        (unless (rczerop sum) (add-term-&-pop a mr sum))
  1602.        (setq m (f1+ m) mr (e+ mr inc))
  1603.        (go b)
  1604.      end
  1605.        (return (pscheck varh (list trunc) (cdr ans)))))
  1606.  
  1607. ;; Computes log(monom), where monom = c x^n. Is extra careful trying to keep
  1608. ;; singular logs going to INF and not generating log(-1)'s unless it is
  1609. ;; necessary to transform a log at MINF to INF.
  1610.  
  1611. (defun pslog-monom (monom)
  1612.   (let* ((gvar (gvar monom))
  1613.      (datum (gvar-data gvar)) var pt logvar c)
  1614.      (if (switch 'MULTIVAR datum) (pslog (ps-lc monom))
  1615.     (setq var (datum-var datum))
  1616.     (if (tvar? var)
  1617.         (if (not (memq (setq pt (exp-pt datum)) '($INF $MINF)))
  1618.         (setq logvar (adjoin-tvar `((%LOG) ,(m- var pt))))
  1619.            ;; At x = inf: log(c (1/x)^n) -> log(c) - n log(x)
  1620.            ;; At x = minf: log(c (-1/x)^n) -> log(c (-1)^n) - n log(x)
  1621.            (setq logvar (psminus (adjoin-tvar `((%LOG) ,var))))
  1622.            (when (eq pt '$MINF) (setq c (rcexpt (rcmone) (ps-le monom)))))
  1623.        (if (eq (caar var) 'MEXPT)
  1624.         (if (equal (caddr var) -1) ;; var must be 1/log(y)
  1625.             ;; Try to keep inf's real. Here we want
  1626.             ;; log(c (1/log(x))^n) -> log(c (-1)^n) - n log(-log(x))
  1627.             (if (equal (tvar-lim (cadr var)) '$MINF)
  1628.             (setq c (rcexpt (rcmone) (ps-le monom))
  1629.                   logvar
  1630.                   (psminus (adjoin-tvar
  1631.                     `((%LOG) ,(m- (cadr var))))))
  1632.                (setq logvar (psminus
  1633.                      (adjoin-tvar `((%LOG) ,(cadr var))))))
  1634.            (if (equal (cadr var) '$%E)
  1635.                (setq logvar (taylor2 (caddr var)))
  1636.               (break "Unhandled gvar in PSLOG-OF-GVAR")))))
  1637.     (psplus (pslog (if c (pstimes c (ps-lc monom)) (ps-lc monom)))
  1638.         (pstimes (ps-le monom) logvar)))))
  1639.  
  1640. ;; Computes log(p), where p is an rcdisrep'd pscoef.
  1641.  
  1642. (defun pslog2 (p) (let ($logarc) (pslog3 p)))
  1643.  
  1644. (defun pslog3 (p)
  1645.    (cond ((atom p)
  1646.       (prep1 (cond ((equal p 1) 0)
  1647.                ((equal p -1) log-1)
  1648.                ((eq p '$%I) log%i)
  1649.                ((eq p '$%E) 1)
  1650.                ((equal p 0)
  1651.             (merror "Log(0) generated while TAYLOR expanding ~M"
  1652.                 last-exp))
  1653.                (T `((%LOG) ,p)))))
  1654.      ((eq (caar p) 'RAT)
  1655.       (prep1 (ifn $TAYLOR_LOGEXPAND `((%LOG) ,p)
  1656.             (m- `((%LOG) ,(cadr p)) `((%LOG) ,(caddr p))))))
  1657.      ((and full-log (not (free p '$%I)))
  1658.       (let ((full-log () )) (pslog3 ($polarform p))))
  1659.      ((eq (caar p) 'MEXPT)
  1660.       ;; Must handle things like x^a, %e^(a*x), etc. which are pscoef's.
  1661.       (pstimes (taylor2 (caddr p)) (pslog (taylor2 (cadr p)))))
  1662.      ((and (eq (caar p) 'MTIMES) $TAYLOR_LOGEXPAND)
  1663.       (do ((l (cddr p) (cdr l))
  1664.            (ans (pslog3 (cadr p)) (psplus ans (pslog3 (car l)))))
  1665.           ((null l) ans)))
  1666.      (T (prep1 `((%LOG) ,p)))))
  1667.  
  1668. ;        (Comment Subtitle Extending Routines)
  1669.  
  1670. (defun getfun-lt (fun)
  1671.    (let ((exp-datum (get (oper-name fun) 'EXP-FORM)))
  1672.     (cond (exp-datum
  1673.            ;; Info not needed yet.
  1674.            ;; (or (atom (car exp-datum))
  1675.            ;;     (setq 0p-funord (copy (cdar exp-datum))))
  1676.            (exp-datum-lt fun exp-datum))
  1677.           ((setq exp-datum (get (oper-name fun) 'SP2))
  1678.            (setq exp-datum (get-lexp (subst (dummy-var) 'SP2VAR exp-datum)
  1679.                      (rcone) ()))
  1680.            ;; Info not needed yet; need to bind lexp-non0 to T when
  1681.            ;; this is used though so n-term will be there.
  1682.            ;; (and (rczerop (le exp-datum))
  1683.            ;;      (setq 0p-funord (le (n-term exp-datum))))
  1684.            (if (psp exp-datum) (ps-lt exp-datum)
  1685.            (term (rczero) exp-datum)))
  1686.           (t (merror "~&~A---Unknown function in getfun-lt~%" fun)))))
  1687.  
  1688. (declare-top (special var))
  1689.  
  1690. (defun getexp-fun (fun var pw)
  1691.   (declare (special var))
  1692.    (let ((exp-datum (copy (get (oper-name fun) 'EXP-FORM))))
  1693.     (cond ((infp pw) (infin-ord-err))
  1694.           ((null exp-datum)
  1695.            (if (null (setq exp-datum
  1696.                    (get-ps-form (if (atom fun) fun (caar fun)))))
  1697.            (merror "~&~A---power series unavailable, internal error~%"
  1698.                fun)
  1699.           (unless (atom fun)
  1700.              (do ((subvals (cdr fun) (cdr subvals))
  1701.               (subs (safe-get (caar fun) 'sp2subs) (cdr subs)))
  1702.              ((or (null subvals) (null subs))
  1703.               (when (or subvals subs)
  1704.                  (merror "~&Incorrect number of subscripts to the ~
  1705.                         DEFTAYLORed function ~A" (caar fun))))
  1706.             (setq exp-datum (MAXIMA-SUBSTITUTE (car subvals) (car subs)
  1707.                             exp-datum))))
  1708.           (ts-formula exp-datum var pw)))
  1709.           ((e> (exp-datum-le fun exp-datum) pw) (pszero var pw))
  1710.           ((setq exp-datum
  1711.              (apply (exp-fun exp-datum)
  1712.                 (if (atom fun) (cons pw (cdr exp-datum))
  1713.                 (cons pw (cons (cdr fun) (cdr exp-datum))))))
  1714.            (cond ((null exp-datum) (pszero var pw))
  1715.              ((psp exp-datum) exp-datum)
  1716.              (t (make-ps var (ncons pw) exp-datum)))))))
  1717.  
  1718. (declare-top(unspecial var))
  1719.               
  1720. (defun EXPEXP-FUNS (pw l sign chng inc)
  1721.        (prog (e lt-l)
  1722.          (setq e (e l) lt-l (setq l (ncons l)))
  1723.         a    (cond ((e> (setq e (e+ e inc)) pw) (return l))
  1724.            (t (add-term-&-pop
  1725.                lt-l
  1726.                e
  1727.                (rctimes (e// sign
  1728.                      (cond ((e= inc (rcone)) e)
  1729.                        ((e* e (e1- e)))))
  1730.                 (cons 1 (cdr (lc lt-l)))))
  1731.               (setq sign (e* sign chng))))
  1732.          (go a)))
  1733.  
  1734. (defun EXPLOG-FUNS (pw l sign chng inc)
  1735.        (prog (e lt-l)
  1736.          (setq e (e l) lt-l (setq l (ncons l)))
  1737.     a    (cond ((e> (setq e (e+ e inc)) pw) (return l))
  1738.            (t (add-term lt-l e (e// sign e))
  1739.               (setq lt-l (n-term lt-l)
  1740.                 sign (e* sign chng))))
  1741.          (go a)))
  1742.  
  1743. (defun EXPTAN-FUNS (pw l chng)
  1744.        (prog (e lt-l sign fact pow)
  1745.          (setq e (e l) lt-l (setq l (ncons l))
  1746.            sign (rcone) fact '(1 . 2) pow '(4 . 1))
  1747.     a    (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l))
  1748.            (t (setq fact (e// fact (e* e (e1+ e)))
  1749.                 pow (e* '(4 . 1) pow)
  1750.                 sign (e* chng sign))
  1751.               (add-term lt-l e (e* (e* sign fact)
  1752.                        (e* (prep1
  1753.                         ($bern (rcdisrep (e1+ e))))
  1754.                            (e* pow (e1- pow)))))
  1755.               (setq lt-l (n-term lt-l))))
  1756.          (go a)))
  1757.  
  1758. (defun EXPCOT-FUNS (pw l sign chng plus)
  1759.        (prog (e lt-l fact pow)
  1760.          (setq e (e l) lt-l (setq l (ncons l))
  1761.            fact (rcone) pow (rcone))
  1762.     a    (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l))
  1763.            (t (setq fact (e// fact (e* e (e1+ e)))
  1764.                 pow (e* '(4 . 1) pow)
  1765.                 sign (e* chng sign))
  1766.               (add-term lt-l e (e* (e* sign fact)
  1767.                        (e* (prep1
  1768.                         ($bern (rcdisrep (e1+ e))))
  1769.                            (e+ pow plus))))
  1770.               (setq lt-l (n-term lt-l))))
  1771.          (go a)))
  1772.  
  1773. (defun EXPSEC-FUNS (pw l chng)
  1774.        (prog (e lt-l sign fact)
  1775.          (setq e (e l) lt-l (setq l (ncons l))
  1776.            sign (rcone)  fact (rcone))
  1777.     a    (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l))
  1778.            (t (setq fact (e// fact (e* e (e1- e)))
  1779.                 sign (e* chng sign))
  1780.               (add-term lt-l e (e* (e* sign fact)
  1781.                        (prep1 ($euler (rcdisrep e)))))
  1782.               (setq lt-l (n-term lt-l))))
  1783.          (go a)))
  1784.  
  1785. (defun EXPASIN-FUNS (pw l chng)
  1786.        (prog (e lt-l sign n d)
  1787.          (setq e (e l) lt-l (setq l (ncons l)) sign 1 n 1 d 1)
  1788.     a    (cond ((e> (setq e (e+ (rctwo) e)) pw) (return l))
  1789.            (t (setq n (times n (car (e- e (rctwo))))
  1790.                 d (times d (car (e1- e)))
  1791.                 sign (times sign chng))
  1792.               (add-term lt-l e  ; need to reduce here ? - check this.
  1793.                 (let ((x (*red (times n sign)
  1794.                           (times d (car e)))))
  1795.                      (if (atom x) x
  1796.                      (cons (cadr x) (caddr x)))))
  1797.               (setq lt-l (n-term lt-l))))
  1798.          (go a)))
  1799.  
  1800. ;;; This is the table of expansion data for known functions.
  1801. ;;; The format of the EXP-FORM property is as follows:
  1802. ;;;    (<name of the expanding routine for the function or
  1803. ;;;      (name . le of n-term) if expansion is of order 0>
  1804. ;;;      <first term in the expansion or the name of a routine which
  1805. ;;;      computes the order when it may depend on parameters (e.g subsripts)>
  1806. ;;;      <data for the expanding routine>)
  1807.  
  1808.  
  1809. (sloop for (fun exp) on
  1810.  '(%EX    ((EXPEXP-FUNS 1 . 1) ((0 . 1) 1 . 1) (1 . 1) (1 . 1) (1 . 1))
  1811.   %SIN   (EXPEXP-FUNS ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1))
  1812.   %COS   ((EXPEXP-FUNS 2 . 1) ((0 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1))
  1813.   %SINH  (EXPEXP-FUNS ((1 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1))
  1814.   %COSH  ((EXPEXP-FUNS 2 . 1) ((0 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1))
  1815.   %LOG   (EXPLOG-FUNS ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (1 . 1))
  1816.   %ATAN  (EXPLOG-FUNS ((1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (2 . 1))
  1817.   %ATANH (EXPLOG-FUNS ((1 . 1) 1 . 1) (1 . 1) (1 . 1) (2 . 1))
  1818.   %COT   (EXPCOT-FUNS ((-1 . 1) 1 . 1) (1 . 1) (-1 . 1) (0 . 1))
  1819.   %CSC   (EXPCOT-FUNS ((-1 . 1) 1 . 1) (-1 . 1) (-1 . 1) (-2 . 1))
  1820.   %CSCH  (EXPCOT-FUNS ((-1 . 1) 1 . 1) (-1 . 1) (1 . 1) (-2 . 1))
  1821.   %COTH  (EXPCOT-FUNS ((-1 . 1) 1 . 1) (1 . 1) (1 . 1) (0 . 1))
  1822.   %TAN   (EXPTAN-FUNS ((1 . 1) 1 . 1) (-1 . 1))
  1823.   %TANH  (EXPTAN-FUNS ((1 . 1) 1 . 1) (1 . 1))
  1824.   %SEC   ((EXPSEC-FUNS 2 . 1) ((0 . 1) 1 . 1) (-1 . 1))
  1825.   %SECH  ((EXPSEC-FUNS 2 . 1) ((0 . 1) 1 . 1) (1 . 1))
  1826.   %ASIN  (EXPASIN-FUNS ((1 . 1) 1 . 1) 1)
  1827.   %ASINH (EXPASIN-FUNS ((1 . 1) 1 . 1) -1)
  1828.   %GAMMA (EXPGAM-FUN ((-1 . 1) 1 . 1))
  1829.   $PSI   (EXPPLYGAM-FUNS plygam-ord))
  1830.   by 'cddr
  1831.   do  (putprop fun exp 'EXP-FORM))
  1832.  
  1833.  
  1834. (defun known-ps (fun) (getl fun '(EXP-FORM SP2)))
  1835.  
  1836. ;            (comment Autoload Properties)
  1837.  
  1838. #+ITS (mautoload '(fasl dsk maxout)
  1839.          (plygam EXPGAM-FUN EXPPLYGAM-FUNS PLYGAM-POLE GAM-CONST
  1840.              PLYGAM-CONST BETA-TRANS))
  1841.  
  1842. ;        (comment Taylor series expansion routines)
  1843.  
  1844. ;;; SRF is only called externally; by RATF and SIMPEXPT.
  1845.  
  1846. (defun srf (x)
  1847.    (let ((exact-poly t) (tlist) (*within-srf?* 'T))
  1848.       (setq x (taylor1 x ()) tlist (mrat-tlist x))
  1849.       ;; Set trunc levels in the local tlist to correspond to the maximum
  1850.       ;; level occuring in any series.
  1851.       (do ((data tlist (cdr data))
  1852.        (truncs (trunc-vector (mrat-ps x) () )))
  1853.       ((null data))
  1854.      (when (and (car truncs) (e> (car truncs) (current-trunc (car data))))
  1855.         (setf (current-trunc (car data)) (car truncs))))
  1856.       x))
  1857.  
  1858. ;;; [var, pt, order, asymp]
  1859.  
  1860. (defmfun $taylor n
  1861.    (if (= n 0) (wna-err '$taylor))
  1862.    (taylor* (arg 1) (listify (f- 1 n))))
  1863.  
  1864. (defun taylor* (arg l)
  1865.    ;; We must bind $MAXTAYORDER to () below because of the problem of constants
  1866.    ;; not retaining their truncation level. This means that when we add a
  1867.    ;; series which has more terms than the user-specified truncation to a
  1868.    ;; constant we must truncate the series with more terms down to the user
  1869.    ;; specified level because, in the worst case, the constant could be a
  1870.    ;; series no better than to the user-specified level. Hence $MAXTAYORDER
  1871.    ;; is essentially useless until the constant problem is fixed. If one
  1872.    ;; decides to not bind $MAXTAYORDER below then the sum routines must
  1873.    ;; be updated to truncate series with more terms than the user-specified
  1874.    ;; level down to that level---taylor(sin(x)^2-cos(x)^2-1,x,0,1) would
  1875.    ;; give x^2+... in this case if the sum routines weren't updated.
  1876.    ;; Also, batch(mquery,160,aljabr) for another truncation bug which crops
  1877.    ;; up when $maxtayorder isn't bound here. Similarly, loadfile(taybad,rl,
  1878.    ;; aljabr) and see tomh's bug note of 4/15/81.
  1879.    (let ((tlist () ) ($maxtayorder () ) (*within-srf?* () )
  1880.      (exact-poly (if l (not $taylor_truncate_polynomials)
  1881.             'User-specified)))
  1882.      (declare(special *within-srf?* ))
  1883.  
  1884.       (parse-tay-args l)
  1885.       (taylor1 arg (ncons tlist))))
  1886.  
  1887. (defun tay-order (n)
  1888.        (let (($float) (modulus))
  1889.       (cond ((eq n '$INF) (ncons (inf)))
  1890.         ((null n) (wna-err '$taylor))
  1891.         ((null (mnump n))
  1892.          (merror "~&~:M---non-numeric expansion order~%" n))
  1893.         (t (ncons (prep1 n))))))
  1894.  
  1895. (defun re-erat (head exp)
  1896.        (taylor1 exp (list (cadddr (cdr head)))))
  1897.  
  1898. (defun parse-tay-args (l)
  1899.    (cond ((null l) )
  1900.      ((or (symbolp (car l)) (not (eq (caaar l) 'MLIST)))
  1901.       (parse-tay-args1 (list (car l) ($ratdisrep (cadr l)) (caddr l)))
  1902.       (parse-tay-args (cdddr l)))
  1903.      ;((or (numberp (car l)) (null (eq (caaar l) 'MLIST)))
  1904.      ; (merror "Variable of expansion not atomic: ~M"
  1905.      ;      (CAR L)))
  1906.      ((do ((l (cddar l) (cdr l)))
  1907.           ((null l) () )
  1908.          (and (or (mnump (car l)) (eq (car l) '$INF))
  1909.           (return 'T)))
  1910.       (parse-tay-args1 (cdar l))
  1911.       (parse-tay-args (cdr l)))
  1912.      (T (parse-tay-args2 (list (car l) (cadr l) (caddr l)))
  1913.         (parse-tay-args (cdddr l)))))
  1914.  
  1915. (defun parse-tay-args1 (l)
  1916.    (if ($listp (car l)) (parse-tay-args2 l)
  1917.       (let ((v (car l))
  1918.         (pt ($ratdisrep (cadr l)))
  1919.         (ord (tay-order (caddr l)))
  1920.         (switches (make-switch-list (cdddr l))))
  1921.      (push (list v ord pt switches) tlist))))
  1922.  
  1923. (defun parse-tay-args2 (l)
  1924.    (let ((label (gensym))
  1925.      (vs (cdar l))
  1926.      (pts (make-long-list (if ($listp (cadr l))
  1927.                   (copy-top-level (cdadr l))
  1928.                  (ncons (ratdisrep (cadr l))))))
  1929.      (ord (caddr l))
  1930.      (switches (make-switch-list (cdddr l)))
  1931.      (lcm 1)
  1932.      (max 0))
  1933.       (if (atom ord)
  1934.       (setq lcm ord max ord ord (make-long-list (ncons ord)))
  1935.      (do ((a vs (cdr a))
  1936.           (l (cdr ord) (cdr l)))
  1937.          ((null a) (setq ord (cdr ord)))
  1938.         (ifn l (merror "Ran out of truncation levels")
  1939.         (setq lcm (lcm lcm (car l))
  1940.               max (max max (car l))))))
  1941.       (push (list label (tay-order max) 0
  1942.           (ncons (list 'MULTIVAR lcm vs)))
  1943.         tlist)
  1944.       (do ((vl vs (cdr vl))
  1945.        (ordl ord (cdr ordl))
  1946.        (ptl pts (cdr ptl)))
  1947.       ((null vl) )
  1948.      (ifn ptl (merror "~&Ran out of matching pts of expansion~%")
  1949.         (push (list (car vl) (tay-order (car ordl)) (car ptl)
  1950.             (cons (list 'MULTI label
  1951.                     (timesk lcm (expta (car ordl) -1)))
  1952.                   switches))
  1953.           tlist)))))
  1954.  
  1955. (defun make-switch-list (l) (mapcar #'(lambda (q) (cons q t)) l))
  1956.  
  1957. (defun make-long-list (q) (nconc q q))
  1958.  
  1959. ;;; This checks to ensure that there isn't more than one set of multi-
  1960. ;;; dependent variables with different orders of expansion, e.g.
  1961. ;;; taylor(u+v+x+y,[u,v],0,[2,3],[x,y],0,[5,7]) is one.
  1962.  
  1963. (defun ratwtsetup (l)
  1964.    (do ((l l (cdr l)) (a) (sw))
  1965.        ((null l) )
  1966.       (when (setq a (switch 'MULTIVAR (car l)))
  1967.      (do ((ll (cadr a) (cdr ll)))
  1968.          ((null ll) )
  1969.         (cond ((equal (cadr (switch 'MULTI (get-datum (car ll)))) 1) )
  1970.           (sw (merror "Can't have two sets of multi dependent ~
  1971.                    variables which require~%different orders ~
  1972.                    of expansion"))
  1973.           ('T (setq sw 'T) (return 'T)))))))
  1974.  
  1975. (defmvar $TAYLOR_ORDER_COEFFICIENTS T
  1976.  "When TRUE, coefficients of taylor series will be ordered canonically.")
  1977.  
  1978. (defun taylor1 (e tlist)
  1979.   (declare (special *within-srf?* ))
  1980.   (setq tlist (tlist-merge (nconc (find-tlists e) tlist)))
  1981.   (prog ($zerobern $simp $algebraic genpairs varlist tvars sing-tvars
  1982.      log-1 log%i ivars key-vars ans full-log last-exp
  1983.      mainvar-datum zerolist taylor_simplifier least_term? tvar-limits)
  1984.     (setq tlist (mapcan #'(lambda (d)
  1985.                 (if (tvar? (datum-var d)) (ncons d)
  1986.                    (push d sing-tvars)
  1987.                    () ))
  1988.                 tlist))
  1989.     (setq $zerobern t $simp t $algebraic t last-exp e least_term? 'T
  1990.           log-1 '((%LOG SIMP) -1) log%i '((%LOG SIMP) $%I)
  1991.           tvars (mapcar 'car tlist) varlist (copy1* tvars))
  1992.          (when $taylor_simplifier
  1993.         (setq taylor_simplifier
  1994.           (if (fboundp $taylor_simplifier) $taylor_simplifier
  1995.              'taylor_simplifier_caller)))
  1996.     ;; Ensure that the expansion points don't depend on the expansion vars.
  1997.     ;; This could cause an infinite loop, e.g. taylor(x,x,x,1).
  1998.     (do ((tl tlist (cdr tl)))
  1999.         ((null tl) )
  2000.        (unless (mfree (exp-pt (car tl)) tvars)
  2001.           (merror
  2002.            "TAYLOR: attempt to expand ~M~&at a point depending on ~M"
  2003.            e (caar tl))))
  2004.     ;; This drastic initialization ensures that ALGEBRAIC, TELLRAT, DISREP,
  2005.     ;; etc. prop's are removed from our gensyms. RATSETUP does not appear
  2006.     ;; to do this correctly, e.g. see ASB's bug of 1/10/83 (MQUERY 17).
  2007.     (mapc #'(lambda (g) (setplist g () )) genvar)
  2008.     (ratsetup varlist genvar)
  2009.     (when (and $taylor_order_coefficients (not *within-srf?*)) (newvar e))
  2010.     (orderpointer varlist)
  2011.     (maplist #'(lambda (q g)
  2012.              (push (cons (car g) (car q)) key-vars)
  2013.              (let ((data (get-datum (car q))))
  2014.             (rplaca q (transform-tvar (car q) data))
  2015.             (push (cons (car g) (car q)) ivars)
  2016.             ;(setf (data-gvar-o data)
  2017.             ;      (cons (car g) (valget (car g))))
  2018.             (rplacd (cdddr data)
  2019.                 (cons (car g) (valget (car g))))))
  2020.          (do ((v varlist (cdr v)))
  2021.              ((eq (car v) (car tvars)) v))
  2022.          (do ((v varlist (cdr v))
  2023.               (g genvar (cdr g)))
  2024.              ((eq (car v) (car tvars)) g)))
  2025.     (setq genpairs (mapcar #'(lambda (y z)
  2026.                    (putprop z y 'DISREP)
  2027.                    (cons y (cons (pget z) 1)))
  2028.                    varlist genvar))
  2029.     (ratwtsetup tlist)
  2030.     (setup-multivar-disrep () )
  2031.     (setq mainvar-datum (car (last tlist)))
  2032. ;    (comment
  2033. ;     (mapc #'(lambda (d)
  2034. ;          (when (and (not (atom (datum-var d)))
  2035. ;                 (memq (caar (datum-var d)) '(MEXPT %LOG)))
  2036. ;             (adjoin-tvar (datum-var d))))
  2037. ;          tlist))
  2038.     (mapc #'(lambda (d) (adjoin-sing-datum d)) sing-tvars)
  2039.     (setq ans (catch 'tay-err (taylor3 e)))
  2040.     (return
  2041.      (if (atom (car ans)) (tay-error (car ans) (cadr ans)) ans))))
  2042.  
  2043. (defun transform-tvar (var data)
  2044.    (if (not (tvar? var)) var
  2045.       (cond ((and (signp e (exp-pt data)) (null (switch '$ASYMP data)))
  2046.          var)    ;Simple case
  2047.         ((memq (exp-pt data) '($INF INFINITY))
  2048.          (m^ var -1))
  2049.         ((eq (exp-pt data) '$MINF)
  2050.          (m- (m^ var -1)))
  2051.         ((let ((temp (m- var (exp-pt data))))
  2052.         (if (switch '$ASYMP data) (m^ temp -1) temp))))))
  2053.  
  2054. (defun taylor_simplifier_caller (exp)
  2055.    (mcall $taylor_simplifier exp))
  2056.  
  2057. (defun taylor_simplify_recurse (ps)
  2058.    (if (pscoefp ps) (taylor2 (funcall taylor_simplifier (rcdisrep ps)))
  2059.       (let ((datum (ps-data ps)) (var () ))
  2060.      ;; We must treat multivars like 1, since they'll reappear again
  2061.      ;; when we call taylor2 on their disrep'd coeff's.
  2062.      (if (switch 'MULTIVAR datum) (setq datum () )
  2063.         (setq var (getdisrep (gvar-o ps)))
  2064.         ;; Don't push pw's < 0, else constants will be truncated
  2065.         (push-pw datum (emax (trunc-lvl ps) (rczero))))
  2066.      (do ((terms (terms ps) (n-term terms))
  2067.           (ans (rczero) (psplus (if (null datum)
  2068.                      (taylor_simplify_recurse (lc terms))
  2069.                        (pstimes (taylor_simplify_recurse
  2070.                          (lc terms))
  2071.      ;; Don't do
  2072.      ;;     (taylor2 (funcall taylor_simplifier
  2073.      ;;              (m^ var (edisrep (le terms)))))
  2074.      ;; causes terms to be lost when inverting. E.g.
  2075.      ;; taylor(log(1+exp(-1/x)),x,0,5) calls psexpt(<exp(1/x)^3>...3,-1)
  2076.      ;; which must return a series good to 3+3(-1-1)=-3 which, when added
  2077.      ;; to other terms will truncate them to degree -3 also.
  2078.                       (if (ezerop (le terms)) (rcone)
  2079.                          (make-ps ps
  2080.                         (ncons
  2081.                          (term (le terms) (rcone)))))))
  2082.  
  2083.                     ans)))
  2084.          ((null terms)
  2085.           (when datum (pop-pw datum))
  2086.           ans)))))
  2087.  
  2088. (defun push-pw (datum pw)
  2089.    (push pw (trunc-stack datum))
  2090.    ;; When changing the truncation on an internal multivar we must also
  2091.    ;; propagate the change to all var's which depend upon it. See WGD's
  2092.    ;; bug report of 9/15/82 which exhibits the necessity of this.
  2093.    (when (setq datum (switch 'MULTIVAR datum))
  2094.       (do ((vars (cadr datum) (cdr vars)))
  2095.       ((null vars) )
  2096.      (push pw (trunc-stack (get-datum (car vars)))))))
  2097.  
  2098. (defun pop-pw (datum)
  2099.    (pop (trunc-stack datum))
  2100.    ;; See the comment above in push-pw; here we must undo the propagation.
  2101.    (when (setq datum (switch 'MULTIVAR datum))
  2102.       (do ((vars (cadr datum) (cdr vars)))
  2103.       ((null vars) )
  2104.      (pop (trunc-stack (get-datum (car vars)))))))
  2105.  
  2106. (defun setup-multivar-disrep (mrat?)
  2107.    (let ((varlist varlist) (genvar genvar) (multivars () ))
  2108.       (when mrat?
  2109.      (setq varlist (mrat-varlist mrat?) genvar (mrat-genvar mrat?)))
  2110.       (mapc #'(lambda (datum)
  2111.          (when (switch 'MULTIVAR datum)
  2112.             (push (car datum) multivars)
  2113.             ;; All genvar's corresponding to internally generated
  2114.             ;; multivars must "disappear" when disrep'd. If this
  2115.             ;; were not done then disrep'ing gx*gt would give x*t
  2116.             ;; which, upon, re-tayloring would give (gx*gt)*gt,
  2117.             ;; where t is the internal multivar for x, and gt, gx
  2118.             ;; are their genvars. An example where this would occur is
  2119.             ;; taylor(sin(x+y),[x],0,f1,[y],0,1).
  2120.             (putprop (int-gvar datum) 1 'DISREP)))
  2121.         (if mrat? (mrat-tlist mrat?) tlist))
  2122.       ;; Here we must substitute 1 for any genvars which depend on multivars.
  2123.       ;; For example, taylor(x^n,[x],0,0) generates a multivar^n.
  2124.       (when multivars
  2125.      (do ((expl varlist (cdr expl))    
  2126.           (gvarl genvar (cdr gvarl)))
  2127.          ((null expl) )
  2128.         (unless (mfree (car expl) multivars)
  2129.            (putprop (car gvarl) 1 'DISREP))))))
  2130.  
  2131. ;; An example showing why this flag is need is given by
  2132. ;; taylor(-exp(exp(-1/x)+2/x),x,0,-1). Without it, tstimes and
  2133. ;; taylor_simplify_recurse end up trunc'ing the -1.
  2134.  
  2135. (defvar trunc-constants? 'T)
  2136.  
  2137. (defun taylor3 (e)
  2138.    (cond ((mbagp e) (cons (car e) (mapcar #'TAYLOR3 (cdr e))))
  2139.      ((and (null tlist) (not (eq exact-poly 'User-specified)))
  2140.       (xcons (prep1 e)
  2141.          (list 'MRAT 'SIMP varlist genvar)))
  2142.      (T (xcons (if (null taylor_simplifier) (taylor2 e)
  2143.               (setq e (taylor2 e))
  2144.               (let ((exact-poly () ) (trunc-constants? () ))
  2145.              (taylor_simplify_recurse e)))
  2146.            (list 'MRAT 'SIMP varlist genvar tlist 'TRUNC)))))
  2147.  
  2148. (defun find-tlists (e) (let (*a*) (findtl1 e) *a*))
  2149.  
  2150. (defun findtl1 (e)
  2151.   (cond ((or (atom e) (mnump e)) )
  2152.     ((eq (caar e) 'MRAT)
  2153.      (when (memq 'TRUNC (car e))
  2154.         (push (mapcar #'(lambda (q) (copy q)) ; copy is a macro!!!
  2155.               (mrat-tlist e))
  2156.           *a*)))
  2157.     (T (mapc #'FINDTL1 (cdr e)))))
  2158.  
  2159. (defun tlist-merge (tlists)
  2160.    (do ((tlists tlists (cdr tlists)) (tlist () ))
  2161.        ((null tlists) (nreverse tlist))
  2162.       (do ((a_tlist (car tlists) (cdr a_tlist)) (temp nil))
  2163.       ((null a_tlist) )
  2164.      (if (null (setq temp (get-datum (datum-var (car a_tlist)) T)))
  2165.          (if (prog2 (setq temp (car a_tlist))
  2166.             (or (tvar? (datum-var temp))
  2167.                 (memq (caar (datum-var temp)) '(MEXPT %LOG))))
  2168.          (push (list (datum-var temp) (trunc-stack temp)
  2169.                  (exp-pt temp) (switches temp))
  2170.                tlist)
  2171.         (merror "Illegal taylor variable ~M" (datum-var temp)))
  2172.         (if $maxtayorder
  2173.         ;; We must take the max truncation level when $maxtayorder
  2174.         ;; is T, cf. JPG's bug of 9/15/82.
  2175.         (when (e> (current-trunc (car a_tlist)) (current-trunc temp))
  2176.            (setf (current-trunc temp) (current-trunc (car a_tlist))))
  2177.            (unless (e> (current-trunc (car a_tlist)) (current-trunc temp))
  2178.           (setf (current-trunc temp) (current-trunc (car a_tlist)))))
  2179.         (unless (alike1 (exp-pt temp) (exp-pt (car a_tlist)))
  2180.            (merror "Cannot combine two expressions expanded at ~
  2181.             different points"))
  2182.         (setf (switches temp)
  2183.           (union* (switches temp) (switches (car a_tlist))))))))
  2184.  
  2185. (defun compattlist (list)
  2186.    (do ((l list (cdr l)))
  2187.        ((null l) T)
  2188.       (or (alike1 (exp-pt (get-datum (datum-var (car l)))) (exp-pt (car l)))
  2189.       (return () ))))
  2190.  
  2191. (defun taylor2  (e)
  2192.  (let ((last-exp e))        ;; lexp-non0 should be bound here when needed
  2193.   (cond ((assolike e tlist) (var-expand e 1 () ))
  2194.     ((or (mnump e) (atom e) (mfree e tvars))
  2195.      (if (e> (rczero) (current-trunc mainvar-datum))
  2196.          (pszero (data-gvar-o mainvar-datum)
  2197.              (current-trunc mainvar-datum))
  2198.         (if (and taylor_simplifier (not (atom e)))
  2199.         (let ((e-simp (prep1 (funcall taylor_simplifier e))))
  2200.            (when (and (rczerop e-simp) (not (memq e-simp zerolist)))
  2201.               (push e zerolist))
  2202.            e-simp)
  2203.            (prep1 e))))
  2204.     ((null (atom (caar e))) (merror "Bad arg TAYLOR2 - internal error"))
  2205.     ((eq (caar e) 'MRAT)
  2206.      (if (and (memq 'TRUNC (car e))
  2207.           (compatvarlist varlist (mrat-varlist e)
  2208.                  genvar (mrat-genvar e))
  2209.           (compattlist (mrat-tlist e)))
  2210.          (pstrunc (cdr e))
  2211.         (let ((exact-poly () )) (re-taylor e))))
  2212.     ((eq (caar e) 'MPLUS) (tsplus (cdr e)))
  2213.     ((eq (caar e) 'MTIMES) (tstimes (cdr e)))
  2214.     ((eq (caar e) 'MEXPT) (tsexpt (cadr e) (caddr e)))
  2215.     ((eq (caar e) '%LOG) (tslog (cadr e)))
  2216.     ((and (or (known-ps (caar e)) (get (caar e) 'TAY-TRANS))
  2217.           (not (memq 'array (cdar e)))
  2218.           (try-expansion (if (cddr e) (cdr e) (cadr e))
  2219.                  (caar e))) )
  2220.     ((and (mqapplyp e)
  2221.           (cond ((get (subfunname e) 'SPEC-TRANS)
  2222.              (funcall (get (subfunname e) 'SPEC-TRANS) e))
  2223.             ((known-ps (subfunname e))
  2224.              (try-expansion (caddr e) (cadr e))))) )
  2225.     ((memq (caar e) '(%SUM %PRODUCT)) (tsprsum (cadr e) (cddr e) (caar e)))
  2226.     ((eq (caar e) '%DERIVATIVE) (tsdiff (cadr e) (cddr e) e))
  2227.     ((or (eq (caar e) '%AT)
  2228.          (do ((l (mapcar 'car tlist) (cdr l)))
  2229.          ((null l) t)
  2230.          (or (free e (car l)) (return ()))))
  2231.      (newsym e))
  2232.     (T (let ((exact-poly () ))    ; Taylor series aren't exact
  2233.           (taylor2 (diff-expand e tlist)))))))
  2234.  
  2235. (defun compatvarlist (a b c d)
  2236.    (cond ((null a) T)
  2237.      ((or (null b) (null c) (null d)) () )
  2238.      ((alike1 (car a) (car b))
  2239.       (if (not (eq (car c) (car d))) ()
  2240.          (compatvarlist (cdr a) (cdr b) (cdr c) (cdr d))))
  2241.      (T (compatvarlist a (cdr b) c (cdr d)))))
  2242.  
  2243.  
  2244.  
  2245. (defun re-taylor (mrat)
  2246.    (let ((old-tlist (mrat-tlist mrat)) (old-varlist (mrat-varlist mrat))
  2247.      (old-genvar (mrat-genvar mrat)) old-ivars)
  2248.      (declare (special old-tlist old-ivars))
  2249.       ;; Put back the old disrpes so rcdisrep's will work correctly.
  2250.       (mapc #'(lambda (g v) (putprop g v 'DISREP)) old-genvar old-varlist)
  2251.       (setup-multivar-disrep mrat)
  2252.       (setq old-ivars (mapcar #'(lambda (g v) (cons g v))
  2253.                   old-genvar old-varlist))
  2254.       (prog1 (re-taylor-recurse (mrat-ps mrat))
  2255.          ;; Restore the correct disreps.
  2256.          (mapc #'(lambda (g v) (putprop g v 'DISREP)) genvar varlist)
  2257.          (setup-multivar-disrep () ))))
  2258.  
  2259. (defun re-taylor-recurse (ps)
  2260.   (declare (special old-tlist old-ivars))
  2261.    (if (not (psp ps)) (taylor2 (rcdisrep ps))
  2262.       (let (var (datum () ))
  2263.      (setq var (cdr (assq (gvar ps) old-ivars)))
  2264.      ;; We must treat multivars like 1, since they'll reappear again
  2265.      ;; when we call taylor2 or var-expand below.
  2266.      (if (switch 'MULTIVAR (zl-ASSOC var old-tlist))
  2267.          (setq var () )
  2268.         (when (setq datum (var-data var))
  2269.            (push-pw datum (trunc-lvl ps))))
  2270.      (prog1
  2271.       (do ((terms (terms ps) (n-term terms))
  2272.            (ans (rczero)
  2273.             (psplus (if (null var) (re-taylor-recurse (lc terms))
  2274.                    (pstimes (re-taylor-recurse (lc terms))
  2275.                     (if datum
  2276.                         (var-expand (car datum)
  2277.                             (edisrep (le terms))
  2278.                             () )
  2279.                        (taylor2
  2280.                         (m^t var (edisrep (le terms)))))))
  2281.                 ans)))
  2282.           ((null terms) ans))
  2283.       (when datum (pop-pw datum))))))
  2284.  
  2285. (defun var-expand (var exp dont-truncate?)
  2286.   (let (($keepfloat) ($float) (modulus))
  2287.      (setq exp (prep1 exp)))        ;; exp must be a rational integer
  2288.   (let ((temp (get-datum var)))
  2289.      (cond ((null temp) (merror "Invalid call to var-expand"))
  2290.        ((switch 'MULTI temp)
  2291.         (psexpt (psplus
  2292.              ;; The reason we call var-expand below instead of taylor2
  2293.              ;; is that we must be sure the call is not truncated to
  2294.              ;; 0 which would cause an error in psexpt if exp < 0.
  2295.              ;; For example, this occured in TAYLOR(X^2/Y,[X,Y],0,2).
  2296.              (pstimes
  2297.               ;; Must ensure that we get back a series truncated
  2298.               ;; to at least what is specified by tlist. This means
  2299.               ;; we'll have to push-pw unless exp>0 since psexpt'n
  2300.               ;; kills (exp-1) terms. The bug that discovered this
  2301.               ;; is taylor(li[2](x+1/2)/x,[x],0,0) missing 2*log(2).
  2302.               (if (not (e> exp (rczero)))
  2303.               (let-pw (get-datum (car (switch 'multi temp)))
  2304.                   (e+ (current-trunc temp) (e- (e1- exp)))
  2305.                  (var-expand (car (switch 'MULTI temp)) 1 'T))
  2306.              (var-expand (car (switch 'MULTI temp)) 1 'T))
  2307.               (cons (list (int-gvar temp) 1 1) 1))
  2308.              (taylor2 (exp-pt temp)))
  2309.             exp))
  2310.        ((signp e (exp-pt temp))
  2311.         (let ((exp>trunc? () ))
  2312.            (if (and (e> exp (current-trunc temp)) (setq exp>trunc? 'T)
  2313.             (not dont-truncate?))
  2314.            (rczero)
  2315.           (make-ps (int-var temp)
  2316.                (ncons (if exact-poly (inf)
  2317.                      (if exp>trunc? exp (current-trunc temp))))
  2318.                (ncons (term (if (switch '$ASYMP temp) (rcminus exp)
  2319.                        exp)
  2320.                     (rcone)))))))
  2321.          ((memq (exp-pt temp) '($INF $MINF $INFINITY))
  2322.           (cond ((switch '$ASYMP temp)
  2323.              (merror
  2324.               "Cannot create an asymptotic expansion at infinity"))
  2325.             ((e> (setq exp (rcminus exp)) (current-trunc temp))
  2326.              (rczero))
  2327.             (t (make-ps (int-var temp)
  2328.                 (ncons (if exact-poly (inf) (current-trunc temp)))
  2329.                 (ncons (term exp
  2330.                          (if (eq (exp-pt temp) '$MINF)
  2331.                          (rcmone)
  2332.                          (rcone))))))))
  2333.          (t (psexpt (psplus
  2334.              (make-ps (int-var temp)
  2335.                   (ncons (if exact-poly (inf) (current-trunc temp)))
  2336.                   (ncons (term (if (switch '$ASYMP temp)
  2337.                            (rcmone)
  2338.                            (rcone))
  2339.                            (rcone))))
  2340.              (taylor2 (exp-pt temp)))
  2341.             exp)))))
  2342.  
  2343. (defun expand (arg func)
  2344.    (or (try-expansion arg func) (exp-pt-err)))
  2345.  
  2346. (defun try-expansion (arg func)
  2347.    (prog (funame funord fun-lc argord psarg arg-trunc temp exact-poly)
  2348.       ;; We bind exact-poly to () since we dont want psexpt retaining
  2349.       ;; higher order terms when subst'ing into series (which aren't exact).
  2350.       ;; Try diff-expanding unknown subsripted functions.
  2351.       (unless (or (atom func) (known-ps (caar func)))
  2352.      (taylor2 (diff-expand `((MQAPPLY) ,func ,arg) tlist)))
  2353.       (when (setq temp (get (setq funame (oper-name func)) 'TAY-TRANS))
  2354.      (return (funcall temp arg func)))
  2355.       (let ((lterm (getfun-lt func)))
  2356.      (setq funord (e lterm) fun-lc (c lterm)))
  2357.     begin-expansion
  2358.       (when (rczerop (or psarg (setq psarg (get-lexp arg (rcone) () ))))
  2359.      (if (e> (rczero) funord)
  2360.          (if (rczerop (setq psarg (get-lexp arg (rcone) 'T)))
  2361.          (tay-depth-err)
  2362.         (go begin-expansion))
  2363.         (return (cond ((setq temp (assq funame TAY-POLE-EXPAND))
  2364.                (funcall (cdr temp) arg psarg func))
  2365.               ((rczerop funord) fun-lc)
  2366.               (T (rczero))))))
  2367.       (when (pscoefp psarg) (setq psarg (taylor2 arg)))
  2368.       (when (pscoefp psarg)
  2369.         (return
  2370.          (cond ((null (mfree (rdis psarg) tvars))
  2371.             (symbolic-expand arg psarg func))
  2372.            ((setq temp (assq funame TAY-POLE-EXPAND))
  2373.             (funcall (cdr temp) arg psarg func))
  2374.            (T (prep1 (simplify
  2375.                   (if (atom func) `((,func) ,(rcdisrep psarg))
  2376.                   `((MQAPPLY) ,func ,(rcdisrep psarg)))))))))
  2377.       (when (e> (rczero) (setq argord (ps-le psarg)))
  2378.      (ifn (memq funame '(%ATAN %ASIN %ASINH %ATANH))
  2379.           (if (e> (rczero) (ps-le* (setq psarg (get-lexp arg (rcone) 'T))))
  2380.           (essen-sing-err)
  2381.          (go begin-expansion))
  2382.         (if (and (eq funame '%ATAN)
  2383.              (eq (asksign-p-or-n (term-disrep (ps-lt psarg) psarg))
  2384.              '$POS))
  2385.         (return (psplus (atrigh arg func) (taylor2 '$%PI)))
  2386.            (return (atrigh arg func)))))
  2387.       (setq temp (t-o-var (gvar psarg)))
  2388.       (when (e> (e* funord argord) temp) (return (rczero)))
  2389.       ;; The following form need not be executed if psarg is really exact.
  2390.       ;; The constant problem does not allow one to determine this now,
  2391.       ;; so we always have to execute this currently.
  2392.       ;; This really should be 
  2393.       ;; (unless (infp (trunc-lvl psarg)) ... )
  2394.       ;; Likewise, the infp checks shouldn't be there; have to assume
  2395.       ;; nothing is exact until constant problem is fixed.
  2396.       (setq arg-trunc (if (and (not (infp (trunc-lvl psarg)))
  2397.                    (e= funord (rcone)))
  2398.               temp
  2399.              (e- temp (e* (e1- funord) argord)))
  2400.         psarg (let-pw (get-datum (get-key-var (gvar psarg)))
  2401.               arg-trunc
  2402.              (if (or (infp (trunc-lvl psarg))
  2403.                  (e> arg-trunc (trunc-lvl psarg)))
  2404.              (taylor2 arg)
  2405.             (pstrunc psarg)))
  2406.         ;; We must recalculate argord since pstrunc may have "picked"
  2407.         ;; a coeff out of a constant monomial; e.g. this occurs in
  2408.         ;; taylor(sin(x+y),x,0,0,y,0,1) where psarg is (Y+...)*X^0+...
  2409.         ;; which truncates to Y+... of order 1.
  2410.         argord (ps-le* psarg))
  2411.     (if (rczerop argord)
  2412.         (cond ((memq funame '(%ATAN %ASIN %ASINH %ATANH))
  2413.            (return (atrigh arg func)))
  2414.           ((setq temp (assq funame const-exp-funs))
  2415.            (return (funcall (cdr temp) arg psarg func)))
  2416.           ((rczerop (ps-le* (setq psarg (get-lexp arg (rcone) 'T))))
  2417.            (return () ))    ; Don't know an addition formula
  2418.           (T (go begin-expansion)))
  2419.        (return
  2420.         (if (mono-term? (terms psarg))
  2421.         (get-series func (current-trunc
  2422.                   (get-datum (get-key-var (gvar psarg))))
  2423.                 (gvar-o psarg) (ps-le psarg) (ps-lc psarg))
  2424.            (setq temp (get-series func (e// temp argord) (gvar-o psarg)
  2425.                       (rcone) (rcone)))
  2426.            (ifn (psp temp) temp (pscsubst1 psarg temp)))))))
  2427.  
  2428. (defun symbolic-expand (arg psarg func) ; should be much stronger
  2429.        arg ;Ignored.
  2430.        (prep1 (simplifya (if (atom func) `((,func) ,(rcdisrep psarg))
  2431.                  `((MQAPPLY) ,func ,(rcdisrep psarg)))
  2432.              () )))
  2433.  
  2434. (defun expand-sing-trig? (arg func)
  2435.    (cond ((memq func *pscirc) (tay-exponentialize arg func))
  2436.      ((memq func *psacirc) (atrigh arg func))
  2437.      (T (essen-sing-err))))
  2438.  
  2439. (defun trig-const (a arg func)
  2440.        (let ((const (ps-lc* arg)) (temp (cdr (assq func trigdisp))))
  2441.         (cond ((and (pscoefp const)
  2442.             (memq func '(%TAN %COT))
  2443.             (multiple-%pi a (srdis const) func)))
  2444.           (temp (funcall temp (setq const (psdisrep const))
  2445.                  (m- a const)))
  2446.           (t (tsexpt `((,(get func 'RECIP)) ,(srdis arg)) -1)))))
  2447.  
  2448. (defun multiple-%pi (a const func)
  2449.   (let (coef)
  2450.     (and (equal ($hipow const '$%PI) 1)
  2451.      ($ratnump (setq coef ($ratcoef const '$%PI 1)))
  2452.      (cond ((numberp coef) (expand (m- a const) func))
  2453.            ((equal (caddr coef) 2)
  2454.         (psminus (expand (m- a const)
  2455.                  (cond ((eq func '%TAN) '%COT)
  2456.                        ((eq func '%COT) '%TAN)
  2457.                        (t (merror "Internal MAXIMA-ERROR in TAYLOR"
  2458.                           ))))))))))
  2459.  
  2460. (setq *pscirc '(%COT %TAN %CSC %SIN %SEC %COS %COTH
  2461.         %TANH %CSCH %SINH %SECH %COSH)
  2462.  
  2463.       *psacirc '(%ACOT %ATAN %ACSC %ASIN %ASEC %ACOS %ACOTH
  2464.                %ATANH %ACSCH %ASINH %ASECH %ACOSH))
  2465.  
  2466. (setq const-exp-funs
  2467.       `((%GAMMA . GAM-CONST) ($PSI . PLYGAM-CONST)
  2468.     . ,(mapcar #'(lambda (q) (cons q 'TRIG-CONST)) *pscirc))
  2469.  
  2470.       trigdisp '((%SIN . PSINA+B) (%COS . PSCOSA+B) (%TAN . PSTANA+B)
  2471.                  (%SINH . PSINHA+B) (%COSH . PSCOSHA+B) (%TANH . PSTANHA+B))
  2472.       
  2473.       tay-pole-expand '((%GAMMA . PLYGAM-POLE) ($PSI . PLYGAM-POLE))
  2474.  
  2475.       tay-const-expand ; !these should be handled by symbolic-expand
  2476.                ; be sure to change this with tay-exponentialize!
  2477.       (append (mapcar #'(lambda (q) (cons q 'TAY-EXPONENTIALIZE)) *PSCIRC)
  2478.           (mapcar #'(lambda (q) (cons q 'TAY-EXPONENTIALIZE)) *PSACIRC)))
  2479.  
  2480. (mapc #'(lambda (q) (putprop q 'ATRIG-TRANS 'TAY-TRANS))
  2481.       '(%ACOS %ACOT %ASEC %ACSC %ACOSH %ACOTH %ASECH %ACSCH))
  2482.  
  2483. (defprop MFACTORIAL FACTORIAL-TRANS TAY-TRANS)
  2484.  
  2485. (defun factorial-trans (arg func)
  2486.   (declare (ignore func))
  2487.   (taylor2 `((%GAMMA) ,(m1+ arg))))
  2488.  
  2489. ;;; Not done properly yet
  2490. ;;;
  2491. ;;; (defprop $BETA BETA-TRANS TAY-TRANS)
  2492.  
  2493. (defun psina+b (a b)
  2494.     (psplus (pstimes (expand a '%SIN) (expand b '%COS))
  2495.             (pstimes (expand a '%COS) (expand b '%SIN))))
  2496.  
  2497. (defun pscosa+b (a b)
  2498.     (psdiff (pstimes (expand a '%COS) (expand b '%COS))
  2499.         (pstimes (expand a '%SIN) (expand b '%SIN))))
  2500.  
  2501. (defun pstana+b (a b)
  2502.     (setq a (expand a '%TAN) b (expand b '%TAN))
  2503.     (pstimes (psplus a b)
  2504.          (psexpt (psdiff (rcone) (pstimes a b))
  2505.              (rcmone))))
  2506.  
  2507. (defun psinha+b (a b)
  2508.     (psplus (pstimes (expand a '%SINH) (expand b '%COSH))
  2509.             (pstimes (expand a '%COSH) (expand b '%SINH))))
  2510.  
  2511. (defun pscosha+b (a b)
  2512.     (psplus (pstimes (expand a '%COSH) (expand b '%COSH))
  2513.         (pstimes (expand a '%SINH) (expand b '%SINH))))
  2514.  
  2515. (defun pstanha+b (a b)
  2516.     (setq a (expand a '%TANH) b (expand b '%TANH))
  2517.     (pstimes (psplus a b)
  2518.          (psexpt (psplus (rcone) (pstimes a b))
  2519.              (rcmone))))
  2520.  
  2521. (defun atrig-trans (arg func)
  2522.        (taylor2
  2523.     (if (memq func '(%ACOS %ACOSH))
  2524.         `((MPLUS)
  2525.           ,HALF%PI
  2526.           ((MTIMES) -1
  2527.             ((,(cdr (assq func
  2528.                       '((%ACOS . %ASIN) (%ACOSH . %ASINH)))))
  2529.              ,arg)))
  2530.         `((,(cdr (assq func '((%ACSC . %ASIN) (%ASEC . %ACOS)
  2531.                   (%ACOT . %ATAN) (%ACSCH . %ASINH)
  2532.                   (%ASECH . %ACOSH) (%ACOTH . %ATANH)))))
  2533.           ,(m^ arg -1)))))
  2534.  
  2535. (defun atrigh (arg func)
  2536.        (let ((full-log t) ($logarc t) (log-1 '((MTIMES) $%I $%PI))
  2537.          (log%i '((MTIMES) ((RAT) 1 2) $%I $%PI)))
  2538.         (taylor2 (simplify `((,func) ,arg)))))
  2539.  
  2540. (defun tay-exponentialize (arg fun) ; !this should be in symbolic-expand!
  2541.        (let (($EXPONENTIALIZE t) ($LOGARC t))
  2542.          (setq arg (meval `((,fun) ,arg))))
  2543.        (taylor2 arg))
  2544.  
  2545. (defun tsplus (l)
  2546.        (do ((l (cdr l) (cdr l))
  2547.         (ans (taylor2 (car l))
  2548.          (psplus ans (taylor2 (car l)))))
  2549.        ((null l) ans)))
  2550.  
  2551. (defun ts-formula (form var pw)
  2552.    (let ((datum (get-datum (get-key-var (car var)))))
  2553.       (let-pw datum pw
  2554.      (taylor2 (subst (get-inverse (car var)) 'SP2VAR form)))))
  2555.  
  2556. (defmacro next-series (l) `(cdadr ,l))
  2557.  
  2558. (defun tstimes-get-pw (l pw)
  2559.    (do ((l l (cdr l)) (vect))
  2560.        ((null l) pw)
  2561.       (setq pw (mapcar #'(lambda (pwq ple) (e+ pwq ple))
  2562.                pw (setq vect (ord-vector (cdar l)))))
  2563.       (rplacd (car l) (cons (cdar l) vect))))
  2564.  
  2565. (defun tstimes-l-mult (a)
  2566.    (do ((l (cdr a) (cdr l)) ($maxtayorder t)
  2567.     (a (car a) (pstimes a (car l))))
  2568.        ((null l) a)))
  2569.  
  2570. (defun mzfree (e l)
  2571.    (do ((l l (cdr l)))
  2572.        ((null l) 'T)
  2573.       (or (zfree e (car l)) (return () ))))
  2574.  
  2575. ;;; The lists posl, negl and  zerl have the following format:
  2576. ;;;
  2577. ;;;   ( (<expression> <expansion> <order vector>) . . . )
  2578.  
  2579. (defun tstimes (l)
  2580.   (*bind* ((funl) (expl) (negl) (zerl) (posl)
  2581.        (pw) (negfl) (temp) (fixl (rcone)))
  2582.     (dolist (fun l)            ;; find the exponentials
  2583.        (if (mexptp fun)
  2584.        (push (if (free (caddr fun) (car tvars)) fun
  2585.             `((MEXPT) $%E ,(m* (caddr fun)
  2586.                        `((%LOG) ,(cadr fun)))))
  2587.          expl)
  2588.       (push fun funl)))
  2589.     (when expl
  2590.        (setq expl (tsexp-comb expl))        ;; simplify exps
  2591.        (setq expl (tsbase-comb expl)))        ;; and again
  2592.     (setq l (nconc expl funl))            ;; now try expanding
  2593.     (let ((trunc-constants? () ))
  2594.        (setq expl (cons 0 (mapcar #'(lambda (exp)
  2595.                        (cons exp (taylor2 exp)))
  2596.                   l))) )
  2597.     ;; EXPL is now of the form (0 ( <form> . <taylor2(form)> ) ...)
  2598.     ;; L points behind the cons considered for destructive updating.
  2599.     (do ((l expl) (tem))
  2600.     ((null (cdr l)) )
  2601.        (cond ((rczerop (cdadr l))
  2602.           ;; Consider taylor((a+1/x)*1/x,x,0,-2). Each factor will be on
  2603.           ;; zerl. Each factor will also appear to have le = 0 since its
  2604.           ;; series is 0, which would fool the get-pw routines below if
  2605.           ;; they tried to handle this case. The easiest fix for now
  2606.           ;; appears to be to always call get-lexp here, killing this:
  2607.           (cond ;((null $maxtayorder)
  2608.             ; (setq zerl (cons (cadr l) zerl))
  2609.             ; (rplacd l (cddr l)))
  2610.             ((rczerop (setq tem (get-lexp (caadr l) (rcone) ())))
  2611.              (return (setq zerl 0)))
  2612.             ('T (setq posl (cons (cons (caadr l) tem) posl))
  2613.             (rplacd l (cddr l)))))
  2614.          ((pscoefp (cdadr l))
  2615.           (cond ((mzfree (caadr l) tvars) ;must be zfree to permit ratfuns
  2616.              (setq fixl (pstimes (cdadr l) fixl))
  2617.              (rplacd l (cddr l)))
  2618.             ((setq zerl (cons (cadr l) zerl))
  2619.              (rplacd l (cddr l)))))
  2620.          ((rczerop (ps-le (cdadr l)))
  2621.           (setq zerl (cons (cadr l) zerl))
  2622.           (rplacd l (cddr l)))
  2623.          ((e> (ps-le (cdadr l)) (rczero))
  2624.           (setq posl (cons (cadr l) posl))
  2625.           (rplacd l (cddr l)))
  2626.          ('T (setq l (cdr l)))))
  2627.     (when (equal zerl 0) (return (rczero)))
  2628.     (setq negl (cdr expl) temp (ord-vector fixl))
  2629.     (mapcar #'(lambda (x) (and (e> (rczero) x) (setq negfl t))) temp)
  2630.     (tstimes-get-pw zerl temp)
  2631.     (setq pw (tstimes-get-pw posl (tstimes-get-pw negl temp)))
  2632.     (if (or negl negfl)
  2633.     (setq posl
  2634.           (mapcar #'(lambda (x)
  2635.                (prog2 (mapcar #'(lambda (datum lel pwl)
  2636.                            (push-pw datum
  2637.                           (e+ (current-trunc datum)
  2638.                               (e- lel pwl))))
  2639.                       tlist (cddr x) pw)
  2640.                   (taylor2 (car x))
  2641.                   (mapcar #'(lambda (datum) (pop-pw datum))
  2642.                       tlist)))
  2643.               (nconc posl zerl negl)))
  2644.        (setq posl (nconc (mapcar 'cadr posl) (mapcar 'cadr zerl)
  2645.              (mapcar 'cadr negl))))
  2646.     (setq posl (tstimes-l-mult posl))
  2647.     (let ((ans (cond ((null fixl) posl)
  2648.              ((null posl) fixl)
  2649.              ('T (pstimes fixl posl)))))
  2650.        (if $maxtayorder ans (pstrunc ans)))))
  2651.  
  2652. ;;; This routine transforms a list of exponentials as follows:
  2653. ;;; 
  2654. ;;;    a^c*b^(n*c) ===> (a*b^n)^c,   where n is free of var.
  2655. ;;;
  2656. ;;; This transformation is only applicable when c is not free of var.
  2657.  
  2658. (defun tsexp-comb (l)    ;; ***** clobbers l ***** 
  2659.    (setq l (cons '* l))
  2660.    (do ((a l) (e))    ;; updated by a rplacd or cdr.
  2661.        ((null (cddr a)) (cdr l))    ;; get rid of the *
  2662.       (rplaca (cddadr a) (setq e ($ratsimp (caddr (cadr a)))))
  2663.       ;; Must delete e^0 lest we divide by the 0 below. RWG's byzero bug
  2664.       ;; of 3/1/78 used to cause this.
  2665.       (if (equal e 0) (rplacd a (cddr a))    
  2666.      (if (mfree (caddr (cadr a)) tvars) (pop a)
  2667.         (do ((b (cddr a) (cdr b)) (n))
  2668.         ((null b) (setq a (cdr a)))
  2669.            (when (mfree (setq n ($ratsimp (m// (caddar b)
  2670.                            (caddr (cadr a)))))
  2671.                 tvars)
  2672.           (rplaca b (list '(MEXPT SIMP)
  2673.                   (m* (cadadr a)
  2674.                       (m^ (cadar b) n))    ;; b^n
  2675.                   (caddr (cadr a))))
  2676.           (rplacd a (cddr a))            ;; delete a^c
  2677.           (return () )))))))
  2678.  
  2679. ;;; This routine transforms a list of exponentials as follows:
  2680. ;;; 
  2681. ;;;    a^b*a^c ===> a^(b+c),
  2682. ;;; 
  2683. ;;; this is only necessary when b and c depend on "var."
  2684.  
  2685. (defun tsbase-comb (l)        ;;; *******clobbers l********
  2686.    (setq l (cons '* l))
  2687.    (do ((a l))            ;;; updated by a rplacd or cdr
  2688.        ((null (cddr a)) (cdr l))
  2689.       (do ((b (cddr a) (cdr b)))
  2690.       ((null b) (pop a))    ;;; did not return early so pop.
  2691.      (when (alike1 (cadar b) (cadadr a))
  2692.         (rplaca b (m^ (cadar b) (m+ (caddar b) (caddr (cadr a)))))
  2693.         (rplacd a (cddr a))
  2694.         (return 'T)))))
  2695.  
  2696. (defun tsexpt (b e)
  2697.    (cond ((and (atom b) (mnump e))
  2698.       (if (get-datum b) (var-expand b e () ) (prep1 (m^ b e))))
  2699.      ((mfree e tvars) (tsexpt1 b e))
  2700.      ((eq '$%E b) (tsexpt-red (list e)))
  2701.      (T (tsexpt-red (list (list '(%LOG) b) e)))))
  2702.  
  2703. (defun tsexpt-red (l)
  2704.    (*bind* ((free) (nfree) (full-log) ($logarc t) (expt) (ps) (e)
  2705.         (log-1 '((MTIMES) $%I $%PI))
  2706.         (log%i '((MTIMES) ((RAT) 1 2) $%I $%PI)))
  2707.        (declare (special e))
  2708.     a  (do ((l l (cdr l)))
  2709.        ((null l) )
  2710.       (cond ((mtimesp (car l)) (setq l (append l (cdar l))))
  2711.         ((mfree (car l) tvars) (push (car l) free))
  2712.         (T (push (car l) nfree))))
  2713.        (cond ((or (cdr nfree) (atom (car nfree))) )
  2714.          ((eq (caaar nfree) '%LOG)
  2715.           (return (tsexpt1 (cadar nfree) (m*l free))))
  2716.          ((memq (caaar nfree) *psacirc)
  2717.           (setq l (ncons (simplifya     ;; simplify after removing simp flag
  2718.                   (cons (ncons (caaar nfree)) (cdar nfree))
  2719.                   () ))
  2720.             nfree (cdr nfree))
  2721.           (go a)))
  2722.        ;; Must have truncs > 0 so that logs in the expt aren't trunc'd.
  2723.        ;; E.g, consider taylor(x^(x-1),x,0,-1).
  2724.        (tlist-mapc d (push-pw d (emax (current-trunc d) (rcone))))
  2725.        (setq ps (taylor2 (setq expt (m*l (append nfree free)))))
  2726.        (tlist-mapc d (pop-pw d))
  2727.        ;; Here we must account for the truncation gain or lossage that
  2728.        ;; is encountered in exp(c*log(x)+y) -> x^c*exp(y).
  2729. ;       (comment
  2730. ;       (when (or (pscoefp ps) (not (rczerop (psterm (terms ps) (rczero)))))
  2731. ;      (do ((data tlist (cdr data))
  2732. ;           (const (srdis (if (pscoefp ps) ps
  2733. ;                (psterm (terms ps) (rczero))))))
  2734. ;          ((null data) (setq ps (taylor2 expt)))
  2735. ;         (setq e (psfind-s (taylor2 (ratcoef const
  2736. ;                         (list '(%LOG) (caar data))))))
  2737. ;         (when (memq (exp-pt (car data)) '($INF $MINF $INFINITY))
  2738. ;        (setq e (e- e)))
  2739. ;         (push-pw (car data)
  2740. ;        (e- (emax (current-trunc (car data)) (rczero)) e)))))
  2741.        (let ((c0 (if (pscoefp ps) ps (psterm (terms ps) (rczero))))
  2742.          e^c0 ord-e^c0)
  2743.       (unless (rczerop c0)
  2744.          (setq ord-e^c0 (ord-vector (setq e^c0 (psexpt-fn c0))))
  2745.          ;; Must emax with 0 so that new singular kernals won't be trunc'd
  2746.          ;; e.g exp(1/x+...) to degree -2 should be exp(-1/x)+...
  2747.          ;; Also try taylor(screwa,x,0,-2).
  2748.          (mapc #'(lambda (d o) (push-pw d (emax (e- (current-trunc d) o)
  2749.                            (rczero))))
  2750.            tlist ord-e^c0)
  2751.          (setq ps (psdiff (taylor2 expt) c0)))
  2752.       (setq ps (psexpt-fn ps))
  2753.       (when e^c0
  2754.          (tlist-mapc d (pop-pw d))
  2755.          (setq ps (pstimes e^c0 ps)))
  2756.       (pstrunc ps))))
  2757.  
  2758. ;; Taylor's b^e, where e is independent of tvars.
  2759.  
  2760. (defun tsexpt1 (b e)
  2761.   (prog (s le pw tb)
  2762.     (setq e (let ((modulus () )) ; Don't mod exponents! See WGM's bug
  2763.            (prog2 (mapcar    ;  of 3/6/83 for an example.
  2764.                #'(lambda (datum)
  2765.                 (push-pw datum
  2766.                    (emax (current-trunc datum) (rczero))))
  2767.                tlist)
  2768.               (taylor2 e)
  2769.               (mapcar #'(lambda (datum) (pop-pw datum)) tlist)))
  2770.           s (psfind-s e)
  2771.           tb (taylor2 b)
  2772.           pw (if (psp tb) (current-trunc (get-datum
  2773.                       (get-key-var (gvar tb))))
  2774.             ;; Constant problem kludge.
  2775.             (if (rczerop tb) (current-trunc (car tlist)) (rczero))))
  2776.     (if (floatp (car s))
  2777.         (setq s (MAXIMA-RATIONALIZE (*quo (car s) (cdr s)))))
  2778.     ;; We must ensure that the lc is non-zero since it will be inverted in
  2779.     ;; psexpt.
  2780.     (setq tb (strip-zeroes tb 'T))
  2781.     (cond ((rczerop tb)
  2782.            (when (or ;; When 1 > s we need more terms since -le*(s-1) > 0.
  2783.                  (e> (rcone) s)
  2784.              (and (e> (rczero) pw) (e> s (rcone))))
  2785.           (setq tb (get-lexp b () 'T)))
  2786.            (setq le (ps-le* tb)))
  2787.           ((psp tb) (setq le (ps-le tb)))
  2788.           (T (return (rcexpt tb e))))
  2789.     (and (e> (e* s le) pw) (null $maxtayorder) (return (rczero)))
  2790.     (setq s (e- pw (e* le (e1- s))))
  2791.     ;(setq le (increment-truncs tb))
  2792.     (return
  2793.      (psexpt
  2794.       (if (e> pw s)
  2795.           (if $maxtayorder tb
  2796.          (pstrunc1 tb (list (cons (gvar tb) s))))
  2797.          ;; because of constants not retaining info, have to
  2798.          ;; just keep the constant here
  2799.          (ifn (psp tb) tb
  2800.         (let-pw (get-datum (get-key-var (gvar tb))) s
  2801.            (strip-zeroes (taylor2 b) 'T))))
  2802.       e))))
  2803.  
  2804. ;;; The method of calculating truncation levels below is incorrect.
  2805. ;;; (i.e. increment-truncs & decrement-truncs, also used above)
  2806. ;;; Examples which exhibit this incorrectness are:
  2807. ;;; taylor(log(sin(y)+x),x,0,2,y,0,1) is missing a y/6*x and -1/6*x^2
  2808. ;;; taylor(log(sin(z)+sin(y)+x),x,0,f1,y,0,3,z,0,5) misses a z^5*y^3 term.
  2809.  
  2810. ;;; TSLOG must find the lowest degree term in the expansion of the
  2811. ;;; log arg, then expand with the orders of all var's in this low term
  2812. ;;; incremented by their order in this low term. Note that this is
  2813. ;;; only necessary for var's with ord > 0, since otherwise we have
  2814. ;;; already expanded to a higher ord than required. Also we must
  2815. ;;; not do this for var's with trunc < 0, since this may incorrectly
  2816. ;;; truncate terms which should end up as logs.
  2817.  
  2818. (defun increment-truncs (ps)
  2819.    (do ((ps ps (ps-lc ps)) (trunc (t-o-var (gvar ps))) (data () ))
  2820.        ((pscoefp ps) data)
  2821.       (when (e> (ps-le ps) (rczero))
  2822.      (push (assq (get-key-var (gvar ps)) tlist) data)
  2823.      (push-pw (car data) (e+ (e* (e+ trunc (rctwo)) (ps-le ps))
  2824.                  (current-trunc (car data))))
  2825.      (setq trunc (e+ trunc (current-trunc (car data))))
  2826.      )))
  2827.  
  2828. (defun decrement-truncs (data)
  2829.    (mapc #'(lambda (data) (pop-pw data)) data))
  2830.  
  2831. (defun tslog (arg)
  2832.   (let ((psarg (taylor2 arg)) datum)
  2833.    (when (rczerop psarg) (setq psarg (get-lexp arg () 'T)))
  2834.    ;; We must ensure that the lc is non-zero since it will be inverted in pslog
  2835.    (setq psarg (strip-zeroes psarg 'T))
  2836.    (do ((ps psarg (ps-lc ps)) (shift (rcone) (e* shift (rctwo))))
  2837.        ((pscoefp ps)
  2838.     (when datum
  2839.        (when (rczerop (setq psarg (taylor2 arg)))
  2840.           (setq psarg (get-lexp arg () 'T)))
  2841.        (mapc #'(lambda (data) (pop-pw data)) datum))
  2842.     (pslog psarg))
  2843.       (push (get-datum (get-key-var (gvar ps))) datum)
  2844.       (if (and (e> (ps-le ps) (rczero))
  2845.            (e> (current-trunc (car datum)) (rczero)))
  2846.       (push-pw (car datum) (e+ (e* shift (ps-le ps))
  2847.                    (current-trunc (car datum))))
  2848.      (pop datum)))))
  2849.  
  2850. ;; When e-start is non-null we start expanding at order e-start, ... , 2^m,
  2851. ;; then 2^m*pow, instead of the normal sequence pow, ... , 2^m*pow
  2852. ;; (where m = $taylordepth, pow = ord of var). This is done because it is
  2853. ;; usually much more efficient for large, non-trivial expansions when we only
  2854. ;; want the lowest order term.
  2855.  
  2856. (defun get-lexp (exp e-start zerocheck?)
  2857.  (if (equal exp 0) (if zerocheck? (tay-depth-err) (rczero))
  2858.  (tlist-mapc d (push-pw d (or e-start (emax (orig-trunc d) (rcone)))))
  2859.  (do ((psexp) (i (f1+ $TAYLORDEPTH) (f1- i)))
  2860.      ((signp e i)
  2861.       (tlist-mapc d (pop-pw d))
  2862.       (if zerocheck? (tay-depth-err)
  2863.      (zero-warn exp)
  2864.      (rczero)))
  2865.      (declare (fixnum i))
  2866.      (cond ((and (rczerop (setq psexp (if zerocheck?
  2867.                       (strip-zeroes (taylor2 exp) 'T)
  2868.                      (taylor2 exp))))
  2869.          (not (memq exp zerolist))) )
  2870.        ;; Info not needed yet.
  2871.        ;; ((and lexp-non0 (rczerop (le (terms psexp)))
  2872.        ;;       (mono-term? (terms psexp))))
  2873.        (T (tlist-mapc d (pop-pw d))
  2874.           (return psexp)))
  2875.      (cond ((and (= i 1) e-start)
  2876.         (setq e-start () i 2)
  2877.         (tlist-mapc d (push-pw d (prog1 (e* (orig-trunc d) (current-trunc d))
  2878.                         (pop-pw d)))))
  2879.        (T (tlist-mapc d (push-pw d (prog1 (e* (rctwo) (current-trunc d))
  2880.                           (pop-pw d)))))))))
  2881.  
  2882. (defun 1p (x) (or (equal x 1) (equal x 1.0)))
  2883.  
  2884. (defun [max-trunc] ()
  2885.    (do ((l tlist (cdr l)) (emax (rczero)))
  2886.        ((null l) (f1+ (// (car emax) (cdr emax))))
  2887.       (when (e> (current-trunc (car l)) emax)
  2888.      (setq emax (orig-trunc (car l))))))
  2889.  
  2890. (defun tsprsum (f l type)
  2891.    (if (mfree f tvars) (newsym f)
  2892.       (let ((li (ncons (car l))) (hi (caddr l)) (lv (ncons (cadr l))) a aa
  2893.         ($maxtayorder () ))        ;; needed to determine when terms are 0
  2894.      (if (and (numberp (car lv)) (numberp hi) (greaterp (car lv) hi))
  2895.          (if (eq type '%SUM) (taylor2 0) (taylor2 1))
  2896.         (if (eq type '%SUM) (setq type () ))
  2897.         (do ((m (f* ([max-trunc]) (^ 2 $taylordepth))) (k 0 (f1+ k))
  2898.          (ans ;(mlet li lv (taylor2 (meval f)))
  2899.               (taylor2 (MAXIMA-SUBSTITUTE (car lv) (car li) f))))
  2900.         ((equal hi (car lv)) ans)
  2901.            (rplaca lv (m1+ (car lv)))
  2902.            ;; A cheap heuristic to catch infinite recursion when 
  2903.            ;; possible, should be improved in the future
  2904.            (if (> k m) (exp-pt-err)
  2905.           (setq a ;(mlet li lv (taylor2 (setq aa (meval f))))
  2906.             (taylor2 (MAXIMA-SUBSTITUTE (car lv) (car li) f))))
  2907.            (if type 
  2908.            (if (and (1p (car a)) (1p (cdr a)) (not (1p aa)))
  2909.                (return ans)
  2910.               (setq ans (pstimes a ans)))
  2911.           (if (and (rczerop a) (not (signp e aa)))
  2912.               (return ans)
  2913.              (setq ans (psplus ans a)))))))))
  2914.  
  2915. (defun tsdiff (e l check)
  2916.     (*bind* ((n) (v) (u))
  2917.           (do ((l l (cddr l)))
  2918.           ((null l))
  2919.           (if (and (atom (car l)) (numberp (cadr l))
  2920.                (assq (car l) tlist))
  2921.               (setq n (cons (cadr l) n) v (cons (car l) v))
  2922.               (setq u (cons (car l) (cons (cadr l) u)))))
  2923.           (or n (return (prep1 check)))
  2924.           (if u (setq e (meval (cons '($DIFF) (cons e l)))))
  2925.           (setq l (mapcar #'(lambda (x) (get-datum x)) v))
  2926.           (mapcar #'(lambda (datum pw)
  2927.               (push-pw datum (e+ (current-trunc datum) (prep1 pw))))
  2928.               l n)
  2929.           (setq e (taylor2 e))
  2930.           (mapc #'(lambda (datum) (pop-pw datum)) l)
  2931.           (do ((vl v (cdr vl))
  2932.            (nl n (cdr nl)))
  2933.           ((null vl ) e)
  2934.           (do ((i 1 (f1+ i)))
  2935.               ((> i (car nl)) )
  2936.               (mapc #'(lambda (a b)
  2937.                 (putprop a (prep1 (sdiff b (car v)))
  2938.                      'DIFF))
  2939.                 genvar varlist)
  2940.               (setq e (psdp e))))))
  2941.  
  2942.  
  2943. (defun no-sing-err (x)            ;; try to catch all singularities
  2944.   (let ((errorsw t))
  2945.     (declare (special errorsw))
  2946.     (let ((ans (catch 'errorsw (eval x))))
  2947.       (if (eq ans t) (unfam-sing-err) ans))))
  2948.  
  2949. (defun check-inf-sing (pt-list) ; don't know behavior of random fun's @ inf
  2950.        (and (or (memq '$inf pt-list) (memq '$minf pt-list))
  2951.         (unfam-sing-err)))
  2952.          
  2953. (defun diff-expand (exp l)    ;l is tlist
  2954.        (check-inf-sing (mapcar (function caddr) l))
  2955.        (ifn l exp
  2956.         (setq exp (diff-expand exp (cdr l)))
  2957.         (do ((deriv (sdiff exp (caar l))
  2958.             (sdiff deriv var))
  2959.          (var (caar l))
  2960.          (coef 1 (times coef (f1+ cnt)))
  2961.          (cnt 1 (f1+ cnt))
  2962.          (pt (exp-pt (car l)))
  2963.          (lim (rcdisrep (current-trunc (car l))))
  2964.          (ans (list
  2965.                (no-sing-err
  2966.             `(meval
  2967.               '(($AT) ,exp
  2968.                   ((MEQUAL) ,(caar l)
  2969.                         ,(exp-pt (car l)))))))
  2970.               (cons
  2971.                `((MTIMES)
  2972.              ((RAT SIMP) 1 ,coef)
  2973.              ,(no-sing-err
  2974.                `(meval '(($AT) ,deriv
  2975.                        ((MEQUAL) ,var ,pt))))
  2976.              ((MEXPT) ,(sub* var pt) ,cnt))
  2977.                ans)))
  2978.         ((or (great cnt lim) (equal deriv 0))
  2979.          (cons '(MPLUS) ans)))))
  2980.  
  2981. ;        (Comment  Subtitle Disreping routines)
  2982.  
  2983. (defun edisrep (e)
  2984.        (if (= (cdr e) 1) (car e) (list '(RAT) (car e) (cdr e))))
  2985.  
  2986. (defun striptimes (a)
  2987.        (if (mtimesp a) (cdr a) (ncons a)))
  2988.  
  2989. (defun srdis (x)
  2990.    (let (($psexpand () )) ; Called only internally, no need to expand.
  2991.       ($ratdisrep
  2992.        (cons (list 'MRAT 'SIMP varlist genvar tlist 'TRUNC)
  2993.          x))))
  2994.  
  2995. (defun srdisrep (r)
  2996.    (let ((varlist (mrat-varlist r)) (genvar (mrat-genvar r)))
  2997.       (mapc #'(lambda (exp genv) (putprop genv exp 'DISREP))
  2998.         varlist genvar)
  2999.       (setup-multivar-disrep r)
  3000.       ;; This used to return 0 if psdisrep returned () but this is wrong
  3001.       ;; since taylor(false,x,0,0) would lose. If psdisrep really wants to
  3002.       ;; return () for 0 then we will probably find out soon.
  3003.       (if (eq $psexpand '$MULTI) (psdisexpand (cdr r))
  3004.      (psdisrep (cdr r)))))
  3005.  
  3006. (defun psdisrep (p)
  3007.    (if (psp p)
  3008.        (psdisrep+ (psdisrep2 (terms p) (getdisrep (gvar-o p)) (trunc-lvl p))
  3009.           (if (or $psexpand (trunc-lvl p)) '(MPLUS TRUNC)
  3010.              '(MPLUS EXACT)))
  3011.       (rcdisrep p)))
  3012.  
  3013. (defun psdisrep^ (n var)
  3014.    ;; If var = () then it is an internal var generated in a multivariate
  3015.    ;; expansion so it shouldn't be displayed. If var = 1 then it probably
  3016.    ;; resulted from the substitution in srdisrep, so it depends on an
  3017.    ;; internal var and likewise shouldn't be displayed.
  3018.    (cond ((or (rczerop n) (null var) (equal var 1)) 1)
  3019.      ((equal n (rcone)) var)
  3020.      ((and ps-bmt-disrep (mexptp var) (equal (caddr var) -1))
  3021.       (psdisrep^ (e- n) (cadr var)))
  3022.      ('T `((MEXPT RATSIMP) ,var ,(edisrep n)))))
  3023.  
  3024. ;;; There used to be a hack below that would print a series consisting
  3025. ;;; of merely one term as exact polynomial (i.e. no trailing "..."'s).
  3026. ;;; This is, of course, wrong but the problem with the fix is that
  3027. ;;; now exact things like taylor(y*x,x,0,f1,y,0,1) will display like
  3028. ;;; (y+...) x+... because of the problem with $MAXTAYORDER being internally
  3029. ;;; bound to ()---which causes exact things to look inexact, such as
  3030. ;;; x and y above. See the comment above taylor* for the $MAXTAYORDER problem.
  3031.  
  3032. (defun psdisrep+ (p plush &aux lowest-degree-term)
  3033.    (if ;; An exact sum of one arg is just that arg.
  3034.        (and (null (cdr p)) (eq (cadr plush) 'EXACT)) (car p)
  3035.       ;; Since the DISPLAY package prints trunc'd sum's arguments
  3036.       ;; from right to left we must put the terms of any constant term
  3037.       ;; in decreasing order. Note that only a constant (wrt to the
  3038.       ;; mainvar) term can be a term which is a sum.
  3039.       (when (mplusp (setq lowest-degree-term (car (last p))))
  3040.      (rplacd lowest-degree-term (nreverse (cdr lowest-degree-term))))
  3041.       (cons plush p)))
  3042.  
  3043. (defun psdisrep* (a b)
  3044.      (cond ((equal a 1) b)
  3045.            ((equal b 1) a)
  3046.            (t (cons '(MTIMES RATSIMP)
  3047.             (nconc (striptimes a) (striptimes b))))))
  3048.  
  3049. (defun psdisrep2 (p var trunc)
  3050.    (if (or $ratexpand $psexpand) (psdisrep2expand p var)
  3051.       (do ((a () (cons (psdisrep* (psdisrep (lc p)) (psdisrep^ (le p) var))
  3052.                a))
  3053.        (p p (cdr p)))
  3054.       ((or (null p) (e> (le p) trunc)) a))))
  3055.  
  3056. (defun psdisrep2expand (p var)
  3057.    (do ((p p (cdr p))
  3058.     (l () (nconc (psdisrep*expand (psdisrep (lc p)) (psdisrep^ (le p) var))
  3059.              l)))
  3060.        ((null p) l)))
  3061.  
  3062. (defun psdisrep*expand (a b)
  3063.   (cond ((equal a 1) (list b))
  3064.     ((equal b 1) (list a))
  3065.     ((null (mplusp a))
  3066.      (list (cons '(MTIMES RATIMES) (nconc (striptimes a) (striptimes b)))))
  3067.     ('T (mapcar #'(lambda (z) (psdisrep* z b))
  3068.             (cdr a)))))
  3069.  
  3070.  
  3071.  
  3072. (defun psdisexpand (p)
  3073.   (let ((ans (ncons ())))
  3074.     (declare (special ans))    ;used in pans-add
  3075.     (psdisexcnt p () (rczero))
  3076.        (setq ans
  3077.          (nreverse
  3078.           (mapcar #'(lambda (x) (ifn (cddr x) (cadr x)
  3079.                      (cons '(MPLUS TRUNC) (cdr x))))
  3080.               (cdr ans))))
  3081.        (ifn (cdr ans) (car ans) (cons '(MPLUS TRUNC) ans))))
  3082.  
  3083. (defun psdisexcnt (p l n)
  3084.   (if (psp p)
  3085.       (do ((var (getdisrep (gvar-o p))) (ll (terms p) (n-term ll)))
  3086.       ((null ll) ())
  3087.     (if (rczerop (le ll)) (psdisexcnt (lc ll) l n)
  3088.         (psdisexcnt (lc ll)
  3089.             (cons (psdisrep^ (le ll) var) l)
  3090.             (e+ (le ll) n))))
  3091.       (psans-add (ifn l (rcdisrep p)
  3092.               (psdisrep* (rcdisrep p)
  3093.                  (ifn (cdr l) (car l)
  3094.                       (cons '(MTIMES TRUNC) l))))
  3095.          n)))
  3096.  
  3097. (defun psans-add (exp n)
  3098.   (declare (special ans))    ;bound in psdisexpand
  3099.   (do ((l ans (cdr l)))
  3100.       ((cond ((null (cdr l)) (rplacd l (ncons (list n exp))))
  3101.          ((e= (caadr l) n) (rplacd (cadr l) (cons exp (cdadr l))))
  3102.          ((e> (caadr l) n) (rplacd l (cons (list n exp) (cdr l))))))))
  3103.  
  3104. (defun srconvert (r)
  3105.   (ifn (atom (caadr (cdddar r))) (cons (car r) (psdisextend (cdr r)))
  3106.        (*bind* ((trunclist (cadr (cdddar r)))
  3107.         (tlist) (gps) (temp) 
  3108.         (vs (caddar r))
  3109.         (gens (cadddr (car r))))
  3110.            (setq gps (mapcar 'cons gens vs))
  3111.            (do ((tl (cdr trunclist) (cddr tl)))
  3112.            ((null tl) (cons (list 'MRAT 'SIMP vs gens tlist 'TRUNC)
  3113.                     (srconvert1 (cdr r))))
  3114.            (setq temp (cdr (assq (car tl) gps)))
  3115.            (cond ((null (memq (car tl) (cdr trunclist))))
  3116.              ((mplusp temp) (merror "FOO"))
  3117.              (t (setq tlist
  3118.                   (cons (list* temp
  3119.                            (tay-order (cadr tl))
  3120.                            0 ()
  3121.                            (cons (car tl)
  3122.                              (symbol-value (car tl))))
  3123.                     tlist))))))))
  3124.  
  3125. (defun srconvert1 (p)
  3126.        (ifn (memq (car p) genvar) p
  3127.         (do ((l (cdr p) (cddr l))
  3128.          (a () (cons (term (prep1 (car l))
  3129.                    (srconvert1 (cadr l)))
  3130.                  a)))
  3131.         ((null l)
  3132.          (make-ps (cons (car p) (symbol-value (car p)))
  3133.               (tay-order (oldget trunclist (car p)))
  3134.               a)))))
  3135.  
  3136. ;        (Comment Subtitle ERROR HANDLING)
  3137.  
  3138. (defun tay-error (msg exp)
  3139.   (if silent-taylor-flag (throw 'taylor-catch ())
  3140.       (IF EXP
  3141.       (merror "TAYLOR~A~%~%~M" MSG EXP)
  3142.       (merror "TAYLOR~A" MSG))))
  3143.  
  3144. (defun exp-pt-err ()
  3145.        (tay-err " unable to expand at a point specified in:"))
  3146.  
  3147. (defun essen-sing-err ()
  3148.        (tay-err " encountered an essential singularity in:"))
  3149.  
  3150. (defun unfam-sing-err ()
  3151.        (tay-err " encountered an unfamiliar singularity in:"))
  3152.  
  3153. (defun infin-ord-err ()
  3154.        (tay-err ": Expansion to infinite order?"))
  3155.  
  3156. (defun tay-depth-err ()
  3157.        (tay-err ": TAYLORDEPTH exceeded while expanding:"))
  3158.  
  3159. ;        (Comment Subtitle TAYLORINFO)
  3160.  
  3161. (defmfun $taylorinfo (x)
  3162.   (ifn (memq 'TRUNC (car x)) ()
  3163.        (cons '(MLIST)
  3164.          (mapcar
  3165.           #'(lambda (q)
  3166.           (nconc
  3167.            `((MLIST) ,(car q) ,(exp-pt q)
  3168.                  ,(let ((tr (current-trunc q)))
  3169.                    (cond ((null tr) '$INF)
  3170.                      ((equal (cdr tr) 1)
  3171.                       (car tr))
  3172.                      (t `((RAT) ,(car tr) ,(cdr tr))))))
  3173.            (mapcar #'(lambda (w)
  3174.                      (list '(MEQUAL) (car w) (cdr w)))
  3175.                (switches q))))
  3176.           (cadddr (cdar x))))))
  3177.  
  3178. ;;; Local Modes:
  3179. ;;; Lisp let-pw Indent:2
  3180. ;;; End:
  3181.